Changeset e7676a in git

Timestamp:

May 7, 2012, 3:55:50 PM (11 years ago)

Author:

Oleksandr Motsak <motsak@…>

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

517530b37f04ee705b7b271243a265748c447cf5dccd6db99cc44c93669e843ddb6e06d2682769b8

Parents:

a034968a1bff5f6e4e9ad126551a2e9a377d9f6060cc1289a79586c1732ba3885997c3cb49bc01d9

Message:

Merge pull request #100 from mmklee/ezgcd_sw

Ezgcd sw

Files:

• Tst/Short/fermat_gcd_1var.res.gz.uu (added)
• Tst/Short/fermat_gcd_1var.stat (added)
• Tst/Short/fermat_gcd_1var.tst (added)
• Tst/Short/fermat_gcd_2var.res.gz.uu (added)
• Tst/Short/fermat_gcd_2var.stat (added)
• Tst/Short/fermat_gcd_2var.tst (added)
• Tst/Short/fermat_gcd_3var.res.gz.uu (added)
• Tst/Short/fermat_gcd_3var.stat (added)
• Tst/Short/fermat_gcd_3var.tst (added)
• Tst/Short/fermat_gcd_4var.res.gz.uu (added)
• Tst/Short/fermat_gcd_4var.stat (added)
• Tst/Short/fermat_gcd_4var.tst (added)
• Tst/Short/fermat_gcd_mod_1var.res.gz.uu (added)
• Tst/Short/fermat_gcd_mod_1var.stat (added)
• Tst/Short/fermat_gcd_mod_1var.tst (added)
• Tst/Short/fermat_gcd_mod_2var.res.gz.uu (added)
• Tst/Short/fermat_gcd_mod_2var.stat (added)
• Tst/Short/fermat_gcd_mod_2var.tst (added)
• Tst/Short/fermat_gcd_mod_3var.res.gz.uu (added)
• Tst/Short/fermat_gcd_mod_3var.stat (added)
• Tst/Short/fermat_gcd_mod_3var.tst (added)
• Tst/Short/fermat_gcd_mod_4var.res.gz.uu (added)
• Tst/Short/fermat_gcd_mod_4var.stat (added)
• Tst/Short/fermat_gcd_mod_4var.tst (added)
• Tst/Short/gcd2test.res.gz.uu (added)
• Tst/Short/gcd2test.stat (added)
• Tst/Short/gcd2test.tst (added)
• Tst/Short/gcd3primtest.res.gz.uu (added)
• Tst/Short/gcd3primtest.stat (added)
• Tst/Short/gcd3primtest.tst (added)
• Tst/Short/gcd3test.res.gz.uu (added)
• Tst/Short/gcd3test.stat (added)
• Tst/Short/gcd3test.tst (added)
• Tst/Short/gcd4test.res.gz.uu (added)
• Tst/Short/gcd4test.stat (added)
• Tst/Short/gcd4test.tst (added)
• Tst/Short/gcd5primtest.res.gz.uu (added)
• Tst/Short/gcd5primtest.stat (added)
• Tst/Short/gcd5primtest.tst (added)
• Tst/Short/gcd5test.res.gz.uu (added)
• Tst/Short/gcd5test.stat (added)
• Tst/Short/gcd5test.tst (added)
• Tst/Short/ok_s.lst (modified) (2 diffs)
• factory/algext.cc (modified) (8 diffs)
• factory/cfNewtonPolygon.cc (modified) (4 diffs)
• factory/cfNewtonPolygon.h (modified) (2 diffs)
• factory/cf_gcd.cc (modified) (11 diffs)
• factory/cf_gcd_smallp.cc (modified) (72 diffs)
• factory/cf_gcd_smallp.h (modified) (3 diffs)
• factory/fac_ezgcd.cc (modified) (3 diffs)

Tst/Short/ok_s.lst

@@ -132 +132 @@
 factorizep_s
 factorizeQa_s
+fermat_gcd_1var
+fermat_gcd_2var
+fermat_gcd_3var
+fermat_gcd_4var
+fermat_gcd_mod_1var
+fermat_gcd_mod_2var
+fermat_gcd_mod_3var
+fermat_gcd_mod_4var
 fetch_s
 fglm1_s
@@ -141 +149 @@
 gaussm_ex
 gcd0_s
+gcd2test
+gcd3primtest
+gcd3test
+gcd4test
+gcd5primtest
+gcd5test
 gcdp_s
 gerhard_1_32003_lp

factory/algext.cc

@@ -1 @@
-#include <factory/factoryconf.h>
+#include "config.h"
 
 #ifndef NOSTREAMIO
@@ -24 +24 @@
 #include "cf_map.h"
 #include "cf_generator.h"
+
+#ifdef HAVE_NTL
+#include "NTLconvert.h"
+#endif
 
 /// compressing two polynomials F and G, M is used for compressing,
@@ -585 +589 @@
 }
 
+#ifdef HAVE_NTL
+static CanonicalForm
+myicontent ( const CanonicalForm & f, const CanonicalForm & c )
+{
+    if (f.isOne() || c.isOne())
+      return 1;
+    if ( f.inBaseDomain() && c.inBaseDomain())
+    {
+      if (c.isZero()) return abs(f);
+      return bgcd( f, c );
+    }
+    else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
+              (f.inCoeffDomain() && c.inBaseDomain()) ||
+              (f.inBaseDomain() && c.inCoeffDomain()))
+    {
+      if (c.isZero()) return abs (f);
+      ZZX NTLf= convertFacCF2NTLZZX (f);
+      ZZX NTLc= convertFacCF2NTLZZX (c);
+      NTLc= GCD (NTLc, NTLf);
+      if (f.inCoeffDomain())
+        return convertNTLZZX2CF(NTLc,f.mvar());
+      else
+        return convertNTLZZX2CF(NTLc,c.mvar());
+    }
+    else
+    {
+        CanonicalForm g = c;
+        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
+            g = myicontent( i.coeff(), g );
+        return g;
+    }
+}
+#endif
+
+CanonicalForm
+myicontent ( const CanonicalForm & f )
+{
+#ifdef HAVE_NTL
+    return myicontent( f, 0 );
+#else
+    return 1;
+#endif
+}
+
 CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
 { // f,g in Q(a)[x1,...,xn]
@@ -612 +660 @@
   f = F * bCommonDen(F);
   g = G * bCommonDen(G);
+  CanonicalForm contf= myicontent (f);
+  CanonicalForm contg= myicontent (g);
+  f /= contf;
+  g /= contg;
+  CanonicalForm gcdcfcg= myicontent (contf, contg);
   Variable a, b;
   if(hasFirstAlgVar(f,a))
@@ -627 +680 @@
       Off( SW_RATIONAL );
       Off( SW_USE_QGCD );
-      tmp = gcd( F, G );
+      tmp = gcdcfcg*gcd( f, g );
       On( SW_USE_QGCD );
       if (off_rational) Off(SW_RATIONAL);
@@ -682 +735 @@
       if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
       setCharacteristic(0);
-      return CanonicalForm(1);
+      return gcdcfcg;
     }
     setCharacteristic(0);
@@ -712 +765 @@
         setReduce(a,true);
         if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
-        return tmp;
+        return tmp*gcdcfcg;
       }
       Off( SW_RATIONAL );
@@ -729 +782 @@
   setReduce(a,true);
   Off( SW_USE_QGCD );
-  D = gcd( f, g );
+  D = gcdcfcg*gcd( f, g );
   On( SW_USE_QGCD );
   if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);

factory/cfNewtonPolygon.cc

@@ -222 +222 @@
 }
 
+int **
+merge (int ** points1, int sizePoints1, int ** points2, int sizePoints2,
+       int& sizeResult)
+{
+  int i, j;
+  sizeResult= sizePoints1+sizePoints2;
+  for (i= 0; i < sizePoints1; i++)
+  {
+    for (j= 0; j < sizePoints2; j++)
+    {
+      if (points1[i][0] != points2[j][0])
+        continue;
+      else
+      {
+        if (points1[i][1] != points2[j][1])
+          continue;
+        else
+        {
+          points2[j][0]= -1;
+          points2[j][1]= -1;
+          sizeResult--;
+        }
+      }
+    }
+  }
+  if (sizeResult == 0)
+    return points1;
+
+  int ** result= new int *[sizeResult];
+  for (i= 0; i < sizeResult; i++)
+    result [i]= new int [2];
+
+  int k= 0;
+  for (i= 0; i < sizePoints1; i++, k++)
+  {
+    result[k][0]= points1[i][0];
+    result[k][1]= points1[i][1];
+  }
+  for (i= 0; i < sizePoints2; i++)
+  {
+    if (points2[i][0] < 0)
+      continue;
+    else
+    {
+      result[k][0]= points2[i][0];
+      result[k][1]= points2[i][1];
+      k++;
+    }
+  }
+  return result;
+}
+
 // assumes a bivariate poly as input
 int ** newtonPolygon (const CanonicalForm& F, int& sizeOfNewtonPoly)
@@ -257 +309 @@
     delete [] points[i];
   delete [] points;
+
+  return result;
+}
+
+// assumes a bivariate polys as input
+int ** newtonPolygon (const CanonicalForm& F, const CanonicalForm& G,
+                      int& sizeOfNewtonPoly)
+{
+  int sizeF= size (F);
+  int ** pointsF= new int* [sizeF];
+  for (int i= 0; i < sizeF; i++)
+    pointsF [i]= new int [2];
+  int j= 0;
+  int * buf;
+  int bufSize;
+  for (CFIterator i= F; i.hasTerms(); i++)
+  {
+    buf= getDegrees (i.coeff(), bufSize);
+    for (int k= 0; k < bufSize; k++, j++)
+    {
+      pointsF [j] [0]= i.exp();
+      pointsF [j] [1]= buf [k];
+    }
+    delete [] buf;
+  }
+
+  int sizeG= size (G);
+  int ** pointsG= new int* [sizeG];
+  for (int i= 0; i < sizeG; i++)
+    pointsG [i]= new int [2];
+  j= 0;
+  for (CFIterator i= G; i.hasTerms(); i++)
+  {
+    buf= getDegrees (i.coeff(), bufSize);
+    for (int k= 0; k < bufSize; k++, j++)
+    {
+      pointsG [j] [0]= i.exp();
+      pointsG [j] [1]= buf [k];
+    }
+    delete [] buf;
+  }
+
+  int sizePoints;
+  int ** points= merge (pointsF, sizeF, pointsG, sizeG, sizePoints);
+
+  int n= polygon (points, sizePoints);
+
+  int ** result= new int* [n];
+  for (int i= 0; i < n; i++)
+  {
+    result [i]= new int [2];
+    result [i] [0]= points [i] [0];
+    result [i] [1]= points [i] [1];
+  }
+
+  sizeOfNewtonPoly= n;
+  for (int i= 0; i < sizeF; i++)
+    delete [] pointsF[i];
+  delete [] pointsF;
+  for (int i= 0; i < sizeG; i++)
+    delete [] pointsG[i];
+  delete [] pointsG;
 
   return result;
@@ -484 +598 @@
 
 CanonicalForm
-compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A)
+compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A, bool computeMA)
 {
   int n;
-  int ** newtonPolyg= newtonPolygon (F, n);
-  convexDense (newtonPolyg, n, M, A);
+  int ** newtonPolyg;
+  if (computeMA)
+  {
+    newtonPolyg= newtonPolygon (F, n);
+    convexDense (newtonPolyg, n, M, A);
+  }
   CanonicalForm result= 0;
   ZZ expX, expY;
@@ -565 +683 @@
   }
 
-  for (int i= 0; i < n; i++)
-    delete [] newtonPolyg [i];
-  delete [] newtonPolyg;
-
-  M= inv (M);
+  if (computeMA)
+  {
+    for (int i= 0; i < n; i++)
+      delete [] newtonPolyg [i];
+    delete [] newtonPolyg;
+    M= inv (M);
+  }
+
   return result;
 }

factory/cfNewtonPolygon.h

@@ -40 +40 @@
                      );
 
+/// compute the convex hull of the support of two bivariate polynomials
+///
+/// @return an array of points in the plane which are the vertices of the convex
+///         hull of the support of F and G
+int ** newtonPolygon (const CanonicalForm& F,///< [in] a bivariate polynomial
+                      const CanonicalForm& G,///< [in] a bivariate polynomial
+                      int& sizeOfNewtonPoly  ///< [in, out] size of the result
+                     );
+
 /// check if @a point is inside a polygon described by points
 ///
@@ -76 +85 @@
 compress (const CanonicalForm& F, ///< [in] compressed, i.e. F.level()==2,
                                   ///< bivariate poly
-          mat_ZZ& inverseM,       ///< [in,out] returns the inverse of M
-          vec_ZZ& A               ///< [in,out] returns translation
+          mat_ZZ& inverseM,       ///< [in,out] returns the inverse of M,
+                                  ///< if computeMA==true, M otherwise
+          vec_ZZ& A,              ///< [in,out] returns translation
+          bool computeMA= true    ///< [in] whether to compute M and A
          );
 

factory/cf_gcd.cc

@@ -757 +757 @@
   else if (!fc_and_gc_Univariate)
   {
-    if (
-    isOn(SW_USE_CHINREM_GCD)
-    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
-    )
-    {
-    #if 0
-      if ( p1 == fc.level() )
-        fc = chinrem_gcd( fc, gc );
-      else
-      {
-        fc = replacevar( fc, Variable(p1), Variable(mp) );
-        gc = replacevar( gc, Variable(p1), Variable(mp) );
-        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
-      }
-    #else
-      fc = chinrem_gcd( fc, gc);
-    #endif
-    }
-    else if ( isOn( SW_USE_EZGCD ) )
-    {
-      if ( pe == 1 )
+    if ( isOn( SW_USE_EZGCD ) )
+    {
+      fc= ezgcd (fc, gc);
+      /*if ( pe == 1 )
         fc = ezgcd( fc, gc );
       else if ( pe > 0 )// no variable at position 1
@@ -791 +774 @@
         gc = swapvar( gc, Variable(pe), Variable(1) );
         fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
+      }*/
+    }
+    else if (
+    isOn(SW_USE_CHINREM_GCD)
+    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
+    )
+    {
+    #if 0
+      if ( p1 == fc.level() )
+        fc = chinrem_gcd( fc, gc );
+      else
+      {
+        fc = replacevar( fc, Variable(p1), Variable(mp) );
+        gc = replacevar( gc, Variable(p1), Variable(mp) );
+        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
       }
+    #else
+      fc = chinrem_gcd( fc, gc);
+    #endif
     }
     else
@@ -1236 +1237 @@
 
   CanonicalForm Dn, test= 0;
+  cl =  gcd (f.lc(),g.lc());
+  CanonicalForm gcdcfcg= gcd (cf, cg);
+  CanonicalForm b= 1;
+  int minCommonDeg= 0;
+  for (i= tmax (f.level(), g.level()); i > 0; i--)
+  {
+    if (degree (f, i) <= 0 || degree (g, i) <= 0)
+      continue;
+    else
+    {
+      minCommonDeg= tmin (degree (g, i), degree (f, i));
+      break;
+    }
+  }
+  if (i == 0)
+    return gcdcfcg;
+  for (; i > 0; i--)
+  {
+    if (degree (f, i) <= 0 || degree (g, i) <= 0)
+      continue;
+    else
+      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
+  }
+  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
+     power (CanonicalForm (2), minCommonDeg);
   bool equal= false;
   i = cf_getNumBigPrimes() - 1;
-  cl =  gcd (f.lc(),g.lc());
-
-  CanonicalForm gcdcfcg= gcd (cf, cg);
+
+  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn;
   //Off (SW_RATIONAL);
   while ( true )
@@ -1246 +1271 @@
     p = cf_getBigPrime( i );
     i--;
-    while ( i >= 0 && mod( cl, p ) == 0 )
+    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
     {
       p = cf_getBigPrime( i );
@@ -1253 +1278 @@
     //printf("try p=%d\n",p);
     setCharacteristic( p );
-    Dp = gcd_poly( mapinto( f ), mapinto( g ) );
+    Dp = GCD_small_p (mapinto (f), mapinto (g), cofp, cogp);
     Dp /=Dp.lc();
+    cofp /= lc (cofp);
+    cogp /= lc (cogp);
     setCharacteristic( 0 );
     dp_deg=totaldegree(Dp);
@@ -1265 +1292 @@
     {
       D = mapinto( Dp );
+      cof= mapinto (cofp);
+      cog= mapinto (cogp);
       d_deg=dp_deg;
       q = p;
@@ -1273 +1302 @@
       {
         chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
+        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
+        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
+        cof= newCof;
+        cog= newCog;
         q = newq;
         D = newD;
@@ -1282 +1315 @@
         q = p;
         D = mapinto( Dp );
+        cof= mapinto (cof);
+        cog= mapinto (cog);
         d_deg=dp_deg;
         test= 0;
@@ -1295 +1330 @@
     {
       Dn= Farey(D,q);
+      cofn= Farey(cof,q);
+      cogn= Farey(cog,q);
       int is_rat= isOn (SW_RATIONAL);
       On (SW_RATIONAL);
       cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
+      cofn *= bCommonDen (cofn);
+      cogn *= bCommonDen (cogn);
       if (!is_rat)
         Off (SW_RATIONAL);
@@ -1306 +1345 @@
         equal= true;
       //Dn /=vcontent(Dn,Variable(1));
-      if (equal && fdivides( Dn, f ) && fdivides( Dn, g ) )
+      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
+          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
       {
         //printf(" -> success\n");
@@ -1325 +1365 @@
 }
 
+

factory/cf_gcd_smallp.cc

@@ -34 +34 @@
 #include "facHensel.h"
 #include "cf_iter.h"
+#include "cfNewtonPolygon.h"
 
 // iinline helper functions:
@@ -49 +50 @@
 TIMING_DEFINE_PRINT(gcd_recursion)
 TIMING_DEFINE_PRINT(newton_interpolation)
+
+bool
+terminationTest (const CanonicalForm& F, const CanonicalForm& G,
+                 const CanonicalForm& coF, const CanonicalForm& coG,
+                 const CanonicalForm& cand)
+{
+  CanonicalForm LCCand= abs (LC (cand));
+  if (LCCand*abs (LC (coF)) == abs (LC (F)))
+  {
+    if (LCCand*abs (LC (coG)) == abs (LC (G)))
+    {
+      if (abs (cand)*abs (coF) == abs (F))
+      {
+        if (abs (cand)*abs (coG) == abs (G))
+          return true;
+      }
+      return false;
+    }
+    return false;
+  }
+  return false;
+}
 
 static const double log2exp= 1.442695041;
@@ -216 +239 @@
 
   return 1;
-}
-
-int
-substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
-{
-  if (F.inCoeffDomain() || G.inCoeffDomain())
-    return 0;
-  Variable x= Variable (1);
-  if (degree (F, x) <= 1 || degree (G, x) <= 1)
-    return 0;
-  CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping is pretty expensive
-  CanonicalForm g= swapvar (G, G.mvar(), x);
-  int sizef= 0;
-  int sizeg= 0;
-  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
-  {
-    if (i.exp() == 1)
-      return 0;
-  }
-  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
-  {
-    if (i.exp() == 1)
-      return 0;
-  }
-  int * expf= new int [sizef];
-  int * expg= new int [sizeg];
-  int j= 0;
-  for (CFIterator i= f; i.hasTerms(); i++, j++)
-  {
-    expf [j]= i.exp();
-  }
-  j= 0;
-  for (CFIterator i= g; i.hasTerms(); i++, j++)
-  {
-    expg [j]= i.exp();
-  }
-
-  int indf= sizef - 1;
-  int indg= sizeg - 1;
-  if (expf[indf] == 0)
-    indf--;
-  if (expg[indg] == 0)
-    indg--;
-
-  if (expg[indg] != expf [indf] || (expg[indg] == 1 && expf[indf] == 1))
-  {
-    delete [] expg;
-    delete [] expf;
-    return 0;
-  }
-  // expf[indg] == expf[indf] after here
-  int result= expg[indg];
-  for (int i= indf - 1; i >= 0; i--)
-  {
-    if (expf [i]%result != 0)
-    {
-      delete [] expg;
-      delete [] expf;
-      return 0;
-    }
-  }
-
-  for (int i= indg - 1; i >= 0; i--)
-  {
-    if (expg [i]%result != 0)
-    {
-      delete [] expg;
-      delete [] expf;
-      return 0;
-    }
-  }
-
-  delete [] expg;
-  delete [] expf;
-  return result;
-}
-
-// substiution
-void
-subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A,
-       CanonicalForm& B, const int d
-      )
-{
-  if (d == 1)
-  {
-    A= F;
-    B= G;
-    return;
-  }
-
-  CanonicalForm C= 0;
-  CanonicalForm D= 0;
-  Variable x= Variable (1);
-  CanonicalForm f= swapvar (F, x, F.mvar());
-  CanonicalForm g= swapvar (G, x, G.mvar());
-  for (CFIterator i= f; i.hasTerms(); i++)
-    C += i.coeff()*power (f.mvar(), i.exp()/ d);
-  for (CFIterator i= g; i.hasTerms(); i++)
-    D += i.coeff()*power (g.mvar(), i.exp()/ d);
-  A= swapvar (C, x, F.mvar());
-  B= swapvar (D, x, G.mvar());
-}
-
-CanonicalForm
-reverseSubst (const CanonicalForm& F, const int d)
-{
-  if (d == 1)
-    return F;
-
-  Variable x= Variable (1);
-  if (degree (F, x) <= 0)
-    return F;
-  CanonicalForm f= swapvar (F, x, F.mvar());
-  CanonicalForm result= 0;
-  for (CFIterator i= f; i.hasTerms(); i++)
-    result += i.coeff()*power (f.mvar(), d*i.exp());
-  return swapvar (result, x, F.mvar());
 }
 
@@ -541 +447 @@
 }
 
+CanonicalForm
+GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
+                  CanonicalForm& coF, CanonicalForm& coG,
+                  Variable & alpha, CFList& l, bool& topLevel);
+
+CanonicalForm
+GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
+                  Variable & alpha, CFList& l, bool& topLevel)
+{
+  CanonicalForm dummy1, dummy2;
+  CanonicalForm result= GCD_Fp_extension (F, G, dummy1, dummy2, alpha, l,
+                                          topLevel);
+  return result;
+}
+
 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
 /// l and topLevel are only used internally, output is monic
@@ -547 +468 @@
 CanonicalForm
 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
+                  CanonicalForm& coF, CanonicalForm& coG,
                   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);
+  if (F.isZero() && degree(G) > 0)
+  {
+    coF= 0;
+    coG= Lc (G);
+    return G/Lc(G);
+  }
+  else if (G.isZero() && degree (F) > 0)
+  {
+    coF= Lc (F);
+    coG= 0;
+    return F/Lc(F);
+  }
+  else if (F.isZero() && G.isZero())
+  {
+    coF= coG= 0;
+    return F.genOne();
+  }
+  if (F.inBaseDomain() || G.inBaseDomain())
+  {
+    coF= F;
+    coG= G;
+    return F.genOne();
+  }
+  if (F.isUnivariate() && fdivides(F, G, coG))
+  {
+    coF= Lc (F);
+    return F/Lc(F);
+  }
+  if (G.isUnivariate() && fdivides(G, F, coF))
+  {
+    coG= Lc (G);
+    return G/Lc(G);
+  }
+  if (F == G)
+  {
+    coF= coG= Lc (F);
+    return F/Lc(F);
+  }
 
   CFMap M,N;
   int best_level= myCompress (A, B, M, N, topLevel);
 
-  if (best_level == 0) return B.genOne();
+  if (best_level == 0)
+  {
+    coF= F;
+    coG= G;
+    return B.genOne();
+  }
 
   A= M(A);
@@ -571 +529 @@
   //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 result= gcd (A, B);
+    coF= N (A/result);
+    coG= N (B/result);
+    return N (result);
+  }
 
   CanonicalForm cA, cB;    // content of A and B
@@ -598 +556 @@
   }
 
+  int sizeNewtonPolyg;
+  int ** newtonPolyg= NULL;
+  mat_ZZ MM;
+  vec_ZZ V;
+  bool compressConvexDense= (ppA.level() == 2 && ppB.level() == 2);
+  if (compressConvexDense)
+  {
+    CanonicalForm bufcA= cA;
+    CanonicalForm bufcB= cB;
+    cA= content (ppA, 1);
+    cB= content (ppB, 1);
+    ppA /= cA;
+    ppB /= cB;
+    gcdcAcB *= gcd (cA, cB);
+    cA *= bufcA;
+    cB *= bufcB;
+    if (ppA.isUnivariate() || ppB.isUnivariate())
+    {
+      if (ppA.level() == ppB.level())
+      {
+        CanonicalForm result= gcd (ppA, ppB);
+        coF= N ((ppA/result)*(cA/gcdcAcB));
+        coG= N ((ppB/result)*(cB/gcdcAcB));
+        return N (result*gcdcAcB);
+      }
+      else
+      {
+        coF= N (ppA*(cA/gcdcAcB));
+        coG= N (ppB*(cB/gcdcAcB));
+        return N (gcdcAcB);
+      }
+    }
+
+    newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
+    convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
+
+    for (int i= 0; i < sizeNewtonPolyg; i++)
+      delete [] newtonPolyg[i];
+    delete [] newtonPolyg;
+
+    ppA= compress (ppA, MM, V, false);
+    ppB= compress (ppB, MM, V, false);
+    MM= inv (MM);
+
+    if (ppA.isUnivariate() && ppB.isUnivariate())
+    {
+      if (ppA.level() == ppB.level())
+      {
+        CanonicalForm result= gcd (ppA, ppB);
+        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
+        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
+        return N (decompress (result, MM, V)*gcdcAcB);
+      }
+      else
+      {
+        coF= N (decompress (ppA, MM, V));
+        coG= N (decompress (ppB, MM, V));
+        return N (gcdcAcB);
+      }
+    }
+  }
+
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
   CanonicalForm gcdlcAlcB;
@@ -604 +624 @@
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB))
+  /*if (fdivides (lcA, lcB))
   {
     if (fdivides (A, B))
@@ -613 +633 @@
     if (fdivides (B, A))
       return G/Lc(G);
-  }
+  }*/
 
   gcdlcAlcB= gcd (lcA, lcB);
@@ -621 +641 @@
   if (d == 0)
   {
-    if (substitute > 1)
-      return N (reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
-  }
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
+    return N(gcdcAcB);
+  }
+
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
   if (d0 < d)
@@ -631 +651 @@
   if (d == 0)
   {
-    if (substitute > 1)
-      return N (reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
+    return N(gcdcAcB);
   }
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
-  CanonicalForm newtonPoly;
+  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
+                coG_m, ppCoF, ppCoG;
 
   newtonPoly= 1;
   m= gcdlcAlcB;
   G_m= 0;
+  coF= 0;
+  coG= 0;
   H= 0;
   bool fail= false;
@@ -650 +672 @@
   CanonicalForm prim_elem, im_prim_elem;
   CFList source, dest;
+  int bound1= degree (ppA, 1);
+  int bound2= degree (ppB, 1);
   do
   {
-    random_element= randomElement (m, V_buf, l, fail);
+    random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
     if (fail)
     {
@@ -667 +691 @@
       {
         m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
-        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
+        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
+        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
+        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
         newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
                              source, dest);
@@ -674 +700 @@
         gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
                             source, dest);
+        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
+        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
         for (CFListIterator i= l; i.hasItem(); i++)
           i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
@@ -693 +721 @@
       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);
+      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
                           source, dest);
@@ -699 +729 @@
       gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
                          source, dest);
+      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
 
       fail= false;
-      random_element= randomElement (m, V_buf, l, fail );
+      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
       DEBOUTLN (cerr, "fail= " << fail);
       CFList list;
       TIMING_START (gcd_recursion);
       G_random_element=
-      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
+      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
+                        coF_random_element, coG_random_element, V_buf,
                         list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion,
@@ -717 +750 @@
       TIMING_START (gcd_recursion);
       G_random_element=
-      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
+      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x),
+                        coF_random_element, coG_random_element, V_buf,
                         list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion,
@@ -732 +766 @@
     if (d0 == 0)
     {
-      if (substitute > 1)
-        return N (reverseSubst (gcdcAcB, substitute));
-      else
-        return N(gcdcAcB);
+      coF= N (ppA*(cA/gcdcAcB));
+      coG= N (ppB*(cB/gcdcAcB));
+      return N(gcdcAcB);
     }
     if (d0 >  d)
@@ -747 +780 @@
     (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
     * G_random_element;
+    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
+                        *coF_random_element;
+    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
+                        *coG_random_element;
 
     if (!G_random_element.inCoeffDomain())
@@ -760 +797 @@
       G_m= 0;
       d= d0;
+      coF_m= 0;
+      coG_m= 0;
     }
 
     TIMING_START (newton_interpolation);
     H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
+    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
+    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
     TIMING_END_AND_PRINT (newton_interpolation,
                           "time for newton interpolation: ");
 
     //termination test
-    if (uni_lcoeff (H) == gcdlcAlcB)
-    {
-      cH= uni_content (H);
+    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
+    {
+      if (gcdlcAlcB.isOne())
+        cH= 1;
+      else
+        cH= uni_content (H);
       ppH= H/cH;
+      CanonicalForm lcppH= gcdlcAlcB/cH;
+      CanonicalForm ccoF= lcA/lcppH;
+      ccoF /= Lc (ccoF);
+      CanonicalForm ccoG= lcB/lcppH;
+      ccoG /= Lc (ccoG);
+      ppCoF= coF/ccoF;
+      ppCoG= coG/ccoG;
       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)
+        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+             terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) || 
+             (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
+        {
+          CFList u, v;
+          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
+          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
+          ppCoF= mapDown (ppCoF, prim_elem, im_prim_elem, alpha, u, v);
+          ppCoF= mapDown (ppCoG, prim_elem, im_prim_elem, alpha, u, v);
+          ppH /= Lc(ppH);
+          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
+          if (compressConvexDense)
           {
-            ppH= reverseSubst (ppH, substitute);
-            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+            ppH= decompress (ppH, MM, V);
+            ppCoF= decompress (ppCoF, MM, V);
+            ppCoG= decompress (ppCoG, MM, V);
           }
+          coF= N ((cA/gcdcAcB)*ppCoF);
+          coG= N ((cB/gcdcAcB)*ppCoG);
           return N(gcdcAcB*ppH);
         }
       }
-      else if (fdivides (ppH, A) && fdivides (ppH, B))
-      {
-        if (substitute > 1)
-        {
-          ppH= reverseSubst (ppH, substitute);
-          gcdcAcB= reverseSubst (gcdcAcB, substitute);
+      else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+                (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+                terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
+                (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
+      {
+        if (compressConvexDense)
+        {
+          ppH= decompress (ppH, MM, V);
+          ppCoF= decompress (ppCoF, MM, V);
+          ppCoG= decompress (ppCoG, MM, V);
         }
+        coF= N ((cA/gcdcAcB)*ppCoF);
+        coG= N ((cB/gcdcAcB)*ppCoG);;
         return N(gcdcAcB*ppH);
       }
@@ -803 +866 @@
 
     G_m= H;
+    coF_m= coF;
+    coG_m= coG;
     newtonPoly= newtonPoly*(x - random_element);
     m= m*(x - random_element);
@@ -848 +913 @@
 }
 
+CanonicalForm
+GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
+        CanonicalForm& coF, CanonicalForm& coG,
+        CFList& l, bool& topLevel);
+
+CanonicalForm
+GCD_GF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
+        bool& topLevel)
+{
+  CanonicalForm dummy1, dummy2;
+  CanonicalForm result= GCD_GF (F, G, dummy1, dummy2, l, topLevel);
+  return result;
+}
+
 /// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
 /// Computer Algebra" by Geddes, Czapor, Labahn
@@ -853 +932 @@
 /// faster field arithmetics, however it might fail if the input is large since
 /// the size of the base field is bounded by 2^16, output is monic
-CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
+CanonicalForm
+GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
+        CanonicalForm& coF, CanonicalForm& coG,
         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);
+  if (F.isZero() && degree(G) > 0)
+  {
+    coF= 0;
+    coG= Lc (G);
+    return G/Lc(G);
+  }
+  else if (G.isZero() && degree (F) > 0)
+  {
+    coF= Lc (F);
+    coG= 0;
+    return F/Lc(F);
+  }
+  else if (F.isZero() && G.isZero())
+  {
+    coF= coG= 0;
+    return F.genOne();
+  }
+  if (F.inBaseDomain() || G.inBaseDomain())
+  {
+    coF= F;
+    coG= G;
+    return F.genOne();
+  }
+  if (F.isUnivariate() && fdivides(F, G, coG))
+  {
+    coF= Lc (F);
+    return F/Lc(F);
+  }
+  if (G.isUnivariate() && fdivides(G, F, coF))
+  {
+    coG= Lc (G);
+    return G/Lc(G);
+  }
+  if (F == G)
+  {
+    coF= coG= Lc (F);
+    return F/Lc(F);
+  }
 
   CFMap M,N;
   int best_level= myCompress (A, B, M, N, topLevel);
 
-  if (best_level == 0) return B.genOne();
+  if (best_level == 0)
+  {
+    coF= F;
+    coG= G;
+    return B.genOne();
+  }
 
   A= M(A);
@@ -878 +995 @@
   //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 result= gcd (A, B);
+    coF= N (A/result);
+    coG= N (B/result);
+    return N (result);
+  }
 
   CanonicalForm cA, cB;    // content of A and B
@@ -905 +1022 @@
   }
 
+  int sizeNewtonPolyg;
+  int ** newtonPolyg= NULL;
+  mat_ZZ MM;
+  vec_ZZ V;
+  bool compressConvexDense= (ppA.level() == 2 && ppB.level() == 2);
+  if (compressConvexDense)
+  {
+    CanonicalForm bufcA= cA;
+    CanonicalForm bufcB= cB;
+    cA= content (ppA, 1);
+    cB= content (ppB, 1);
+    ppA /= cA;
+    ppB /= cB;
+    gcdcAcB *= gcd (cA, cB);
+    cA *= bufcA;
+    cB *= bufcB;
+    if (ppA.isUnivariate() || ppB.isUnivariate())
+    {
+      if (ppA.level() == ppB.level())
+      {
+        CanonicalForm result= gcd (ppA, ppB);
+        coF= N ((ppA/result)*(cA/gcdcAcB));
+        coG= N ((ppB/result)*(cB/gcdcAcB));
+        return N (result*gcdcAcB);
+      }
+      else
+      {
+        coF= N (ppA*(cA/gcdcAcB));
+        coG= N (ppB*(cB/gcdcAcB));
+        return N (gcdcAcB);
+      }
+    }
+
+    newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
+    convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
+
+    for (int i= 0; i < sizeNewtonPolyg; i++)
+      delete [] newtonPolyg[i];
+    delete [] newtonPolyg;
+
+    ppA= compress (ppA, MM, V, false);
+    ppB= compress (ppB, MM, V, false);
+    MM= inv (MM);
+
+    if (ppA.isUnivariate() && ppB.isUnivariate())
+    {
+      if (ppA.level() == ppB.level())
+      {
+        CanonicalForm result= gcd (ppA, ppB);
+        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
+        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
+        return N (decompress (result, MM, V)*gcdcAcB);
+      }
+      else
+      {
+        coF= N (decompress (ppA, MM, V));
+        coG= N (decompress (ppB, MM, V));
+        return N (gcdcAcB);
+      }
+    }
+  }
+
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
   CanonicalForm gcdlcAlcB;
@@ -911 +1090 @@
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB))
-  {
-    if (fdivides (A, B))
+  /*if (fdivides (lcA, lcB))
+  {
+    if (fdivides (ppA, ppB, coG))
+    {
+      coF= 1;
+      if (compressConvexDense)
+        coG= decompress (coG, MM, V);
+      coG= N (coG*(cB/gcdcAcB));
       return F;
+    }
   }
   if (fdivides (lcB, lcA))
   {
-    if (fdivides (B, A))
+    if (fdivides (ppB, ppA, coF))
+    {
+      coG= 1;
+      if (compressConvexDense)
+        coF= decompress (coF, MM, V);
+      coF= N (coF*(cA/gcdcAcB));
       return G;
-  }
+    }
+  }*/
 
   gcdlcAlcB= gcd (lcA, lcB);
@@ -927 +1118 @@
   if (d == 0)
   {
-    if (substitute > 1)
-      return N (reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
+    return N(gcdcAcB);
   }
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
@@ -937 +1127 @@
   if (d == 0)
   {
-    if (substitute > 1)
-      return N (reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
+    return N(gcdcAcB);
   }
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
-  CanonicalForm newtonPoly;
+  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
+                coG_m, ppCoF, ppCoG;
 
   newtonPoly= 1;
   m= gcdlcAlcB;
   G_m= 0;
+  coF= 0;
+  coG= 0;
   H= 0;
   bool fail= false;
-  topLevel= false;
+  //topLevel= false;
   bool inextension= false;
   int p=-1;
@@ -958 +1150 @@
   int expon;
   char gf_name_buf= gf_name;
+  int bound1= degree (ppA, 1);
+  int bound2= degree (ppB, 1);
   do
   {
-    random_element= GFRandomElement (m, l, fail);
+    random_element= GFRandomElement (m*lcA*lcB, l, fail);
     if (fail)
     {
@@ -981 +1175 @@
       G_m= GFMapUp (G_m, kk);
       newtonPoly= GFMapUp (newtonPoly, kk);
+      coF_m= GFMapUp (coF_m, kk);
+      coG_m= GFMapUp (coG_m, kk);
       ppA= GFMapUp (ppA, kk);
       ppB= GFMapUp (ppB, kk);
       gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
-      random_element= GFRandomElement (m, l, fail);
+      lcA= GFMapUp (lcA, kk);
+      lcB= GFMapUp (lcB, kk);
+      random_element= GFRandomElement (m*lcA*lcB, l, fail);
       DEBOUTLN (cerr, "fail= " << fail);
       CFList list;
       TIMING_START (gcd_recursion);
       G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
+                                coF_random_element, coG_random_element,
                                 list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion,
@@ -999 +1198 @@
       TIMING_START (gcd_recursion);
       G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
+                                coF_random_element, coG_random_element,
                                 list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion,
@@ -1014 +1214 @@
     {
       if (inextension)
-      {
         setCharacteristic (p, k, gf_name_buf);
-        {
-          if (substitute > 1)
-            return N (reverseSubst (gcdcAcB, substitute));
-          else
-            return N(gcdcAcB);
-        }
-      }
-      else
-      {
-        if (substitute > 1)
-          return N (reverseSubst (gcdcAcB, substitute));
-        else
-          return N(gcdcAcB);
-      }
+      coF= N (ppA*(cA/gcdcAcB));
+      coG= N (ppB*(cB/gcdcAcB));
+      return N(gcdcAcB);
     }
     if (d0 >  d)
@@ -1041 +1229 @@
     (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
      G_random_element;
+
+    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
+                        *coF_random_element;
+    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
+                        *coG_random_element;
+
     if (!G_random_element.inCoeffDomain())
       d0= totaldegree (G_random_element, Variable(2),
@@ -1053 +1247 @@
       G_m= 0;
       d= d0;
+      coF_m= 0;
+      coG_m= 0;
     }
 
     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: ");
+    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
+    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
+    TIMING_END_AND_PRINT (newton_interpolation,
+                          "time for newton interpolation: ");
 
     //termination test
-    if (uni_lcoeff (H) == gcdlcAlcB)
-    {
-      cH= uni_content (H);
+    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
+    {
+      if (gcdlcAlcB.isOne())
+        cH= 1;
+      else
+        cH= uni_content (H);
       ppH= H/cH;
+      CanonicalForm lcppH= gcdlcAlcB/cH;
+      CanonicalForm ccoF= lcA/lcppH;
+      ccoF /= Lc (ccoF);
+      CanonicalForm ccoG= lcB/lcppH;
+      ccoG /= Lc (ccoG);
+      ppCoF= coF/ccoF;
+      ppCoG= coG/ccoG;
       if (inextension)
       {
-        if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k)))
+        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+             terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
+             (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
         {
           DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
           ppH= GFMapDown (ppH, k);
+          ppCoF= GFMapDown (ppCoF, k);
+          ppCoG= GFMapDown (ppCoG, k);
           DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
-          if (substitute > 1)
+          if (compressConvexDense)
           {
-            ppH= reverseSubst (ppH, substitute);
-            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+            ppH= decompress (ppH, MM, V);
+            ppCoF= decompress (ppCoF, MM, V);
+            ppCoG= decompress (ppCoG, MM, V);
           }
+          coF= N ((cA/gcdcAcB)*ppCoF);
+          coG= N ((cB/gcdcAcB)*ppCoG);
           setCharacteristic (p, k, gf_name_buf);
           return N(gcdcAcB*ppH);
@@ -1082 +1299 @@
       else
       {
-        if (fdivides (ppH, A) && fdivides (ppH, B))
-        {
-          if (substitute > 1)
+      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+           terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
+           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
+        {
+          if (compressConvexDense)
           {
-            ppH= reverseSubst (ppH, substitute);
-            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+            ppH= decompress (ppH, MM, V);
+            ppCoF= decompress (ppCoF, MM, V);
+            ppCoG= decompress (ppCoG, MM, V);
           }
+          coF= N ((cA/gcdcAcB)*ppCoF);
+          coG= N ((cB/gcdcAcB)*ppCoG);
           return N(gcdcAcB*ppH);
         }
@@ -1095 +1318 @@
 
     G_m= H;
+    coF_m= coF;
+    coG_m= coG;
     newtonPoly= newtonPoly*(x - random_element);
     m= m*(x - random_element);
@@ -1152 +1377 @@
 }
 
-CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
-                           bool& topLevel, CFList& l)
+CanonicalForm
+GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
+             CanonicalForm& coF, CanonicalForm& coG,
+             bool& topLevel, CFList& l);
+
+CanonicalForm
+GCD_small_p (const CanonicalForm& F, const CanonicalForm& G,
+             bool& topLevel, CFList& l)
+{
+  CanonicalForm dummy1, dummy2;
+  CanonicalForm result= GCD_small_p (F, G, dummy1, dummy2, topLevel, l);
+  return result;
+}
+
+CanonicalForm
+GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
+             CanonicalForm& coF, CanonicalForm& coG,
+             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);
-
+  if (F.isZero() && degree(G) > 0)
+  {
+    coF= 0;
+    coG= Lc (G);
+    return G/Lc(G);
+  }
+  else if (G.isZero() && degree (F) > 0)
+  {
+    coF= Lc (F);
+    coG= 0;
+    return F/Lc(F);
+  }
+  else if (F.isZero() && G.isZero())
+  {
+    coF= coG= 0;
+    return F.genOne();
+  }
+  if (F.inBaseDomain() || G.inBaseDomain())
+  {
+    coF= F;
+    coG= G;
+    return F.genOne();
+  }
+  if (F.isUnivariate() && fdivides(F, G, coG))
+  {
+    coF= Lc (F);
+    return F/Lc(F);
+  }
+  if (G.isUnivariate() && fdivides(G, F, coF))
+  {
+    coG= Lc (G);
+    return G/Lc(G);
+  }
+  if (F == G)
+  {
+    coF= coG= Lc (F);
+    return F/Lc(F);
+  }
   CFMap M,N;
   int best_level= myCompress (A, B, M, N, topLevel);
 
-  if (best_level == 0) return B.genOne();
+  if (best_level == 0)
+  {
+    coF= F;
+    coG= G;
+    return B.genOne();
+  }
 
   A= M(A);
@@ -1177 +1453 @@
   //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 result= gcd (A, B);
+    coF= N (A/result);
+    coG= N (B/result);
+    return N (result);
+  }
 
   CanonicalForm cA, cB;    // content of A and B
@@ -1204 +1480 @@
   }
 
+  int sizeNewtonPolyg;
+  int ** newtonPolyg= NULL;
+  mat_ZZ MM;
+  vec_ZZ V;
+  bool compressConvexDense= (ppA.level() == 2 && ppB.level() == 2);
+  if (compressConvexDense)
+  {
+    CanonicalForm bufcA= cA;
+    CanonicalForm bufcB= cB;
+    cA= content (ppA, 1);
+    cB= content (ppB, 1);
+    ppA /= cA;
+    ppB /= cB;
+    gcdcAcB *= gcd (cA, cB);
+    cA *= bufcA;
+    cB *= bufcB;
+    if (ppA.isUnivariate() || ppB.isUnivariate())
+    {
+      if (ppA.level() == ppB.level())
+      {
+        CanonicalForm result= gcd (ppA, ppB);
+        coF= N ((ppA/result)*(cA/gcdcAcB));
+        coG= N ((ppB/result)*(cB/gcdcAcB));
+        return N (result*gcdcAcB);
+      }
+      else
+      {
+        coF= N (ppA*(cA/gcdcAcB));
+        coG= N (ppB*(cB/gcdcAcB));
+        return N (gcdcAcB);
+      }
+    }
+
+    newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
+    convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
+
+    for (int i= 0; i < sizeNewtonPolyg; i++)
+      delete [] newtonPolyg[i];
+    delete [] newtonPolyg;
+
+    ppA= compress (ppA, MM, V, false);
+    ppB= compress (ppB, MM, V, false);
+    MM= inv (MM);
+
+    if (ppA.isUnivariate() && ppB.isUnivariate())
+    {
+      if (ppA.level() == ppB.level())
+      {
+        CanonicalForm result= gcd (ppA, ppB);
+        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
+        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
+        return N (decompress (result, MM, V)*gcdcAcB);
+      }
+      else
+      {
+        coF= N (decompress (ppA, MM, V));
+        coG= N (decompress (ppB, MM, V));
+        return N (gcdcAcB);
+      }
+    }
+  }
+
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
   CanonicalForm gcdlcAlcB;
@@ -1209 +1547 @@
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB))
+  /*if (fdivides (lcA, lcB))
   {
     if (fdivides (A, B))
@@ -1218 +1556 @@
     if (fdivides (B, A))
       return G/Lc(G);
-  }
+  }*/
 
   gcdlcAlcB= gcd (lcA, lcB);
@@ -1227 +1565 @@
   if (d == 0)
   {
-    if (substitute > 1)
-      return N (reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
-  }
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
+    return N(gcdcAcB);
+  }
+
   d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
 
@@ -1239 +1577 @@
   if (d == 0)
   {
-    if (substitute > 1)
-      return N (reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
+    return N(gcdcAcB);
   }
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
-  CanonicalForm newtonPoly= 1;
+  CanonicalForm newtonPoly, coF_random_element, coG_random_element,
+                coF_m, coG_m, ppCoF, ppCoG;
+
+  newtonPoly= 1;
   m= gcdlcAlcB;
   H= 0;
+  coF= 0;
+  coG= 0;
   G_m= 0;
   Variable alpha, V_buf;
@@ -1254 +1596 @@
   bool inextension= false;
   bool inextensionextension= false;
-  topLevel= false;
+  bool topLevel2= false;
   CFList source, dest;
+  int bound1= degree (ppA, 1);
+  int bound2= degree (ppB, 1);
   do
   {
     if (inextension)
-      random_element= randomElement (m, V_buf, l, fail);
+      random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
     else
-      random_element= FpRandomElement (m, l, fail);
+      random_element= FpRandomElement (m*lcA*lcB, l, fail);
 
     if (!fail && !inextension)
@@ -1268 +1612 @@
       TIMING_START (gcd_recursion);
       G_random_element=
-      GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel,
-      list);
+      GCD_small_p (ppA (random_element,x), ppB (random_element,x),
+                   coF_random_element, coG_random_element, topLevel2,
+                   list);
       TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
@@ -1279 +1624 @@
       TIMING_START (gcd_recursion);
       G_random_element=
-      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
-                        list, topLevel);
+      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
+                        coF_random_element, coG_random_element, alpha,
+                        list, topLevel2);
       TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
@@ -1297 +1643 @@
         inextension= true;
         fail= false;
-        random_element= randomElement (m, alpha, l, fail);
+        random_element= randomElement (m*lcA*lcB, alpha, l, fail);
         deg++;
       } while (fail);
@@ -1304 +1650 @@
       TIMING_START (gcd_recursion);
       G_random_element=
-      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
-                        list, topLevel);
+      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
+                        coF_random_element, coG_random_element, alpha,
+                        list, topLevel2);
       TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
@@ -1325 +1672 @@
       {
         m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
-        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
+        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
+        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
+        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
         newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
                              source, dest);
@@ -1332 +1681 @@
         gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
                             dest);
+        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
+        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
         for (CFListIterator i= l; i.hasItem(); i++)
           i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
@@ -1352 +1703 @@
       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);
+      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
                           source, dest);
@@ -1358 +1711 @@
       gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
                          source, dest);
+      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       fail= false;
-      random_element= randomElement (m, V_buf, l, fail );
+      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
       DEBOUTLN (cerr, "fail= " << fail);
       CFList list;
       TIMING_START (gcd_recursion);
       G_random_element=
-      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
-                        list, topLevel);
+      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
+                        coF_random_element, coG_random_element, V_buf,
+                        list, topLevel2);
       TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
@@ -1379 +1735 @@
     if (d0 == 0)
     {
-      if (substitute > 1)
-        return N (reverseSubst (gcdcAcB, substitute));
-      else
-        return N(gcdcAcB);
-    }
+      coF= N (ppA*(cA/gcdcAcB));
+      coG= N (ppB*(cB/gcdcAcB));
+      return N(gcdcAcB);
+    }
+
     if (d0 >  d)
     {
@@ -1394 +1750 @@
                        *G_random_element;
 
+    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
+                        *coF_random_element;
+    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
+                        *coG_random_element;
 
     if (!G_random_element.inCoeffDomain())
@@ -1407 +1767 @@
       G_m= 0;
       d= d0;
+      coF_m= 0;
+      coG_m= 0;
     }
 
     TIMING_START (newton_interpolation);
     H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
+    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
+    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
     TIMING_END_AND_PRINT (newton_interpolation,
                           "time for newton_interpolation: ");
 
     //termination test
-    if (uni_lcoeff (H) == gcdlcAlcB)
-    {
-      cH= uni_content (H);
+    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
+    {
+      if (gcdlcAlcB.isOne())
+        cH= 1;
+      else
+        cH= uni_content (H);
       ppH= H/cH;
       ppH /= Lc (ppH);
+      CanonicalForm lcppH= gcdlcAlcB/cH;
+      CanonicalForm ccoF= lcppH/Lc (lcppH);
+      CanonicalForm ccoG= lcppH/Lc (lcppH);
+      ppCoF= coF/ccoF;
+      ppCoG= coG/ccoG;
       DEBOUTLN (cerr, "ppH= " << ppH);
-      if (fdivides (ppH, A) && fdivides (ppH, B))
-      {
-        if (substitute > 1)
-        {
-          ppH= reverseSubst (ppH, substitute);
-          gcdcAcB= reverseSubst (gcdcAcB, substitute);
+      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+           terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
+           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
+      {
+        if (compressConvexDense)
+        {
+          ppH= decompress (ppH, MM, V);
+          ppCoF= decompress (ppCoF, MM, V);
+          ppCoG= decompress (ppCoG, MM, V);
         }
+        coF= N ((cA/gcdcAcB)*ppCoF);
+        coG= N ((cB/gcdcAcB)*ppCoG);
         return N(gcdcAcB*ppH);
       }
@@ -1433 +1811 @@
 
     G_m= H;
+    coF_m= coF;
+    coG_m= coG;
     newtonPoly= newtonPoly*(x - random_element);
     m= m*(x - random_element);
@@ -2855 +3235 @@
     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
@@ -2901 +3276 @@
 
   if (d == 0)
-  {
-    if (substitute > 1)
-      return N(reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
-  }
+    return N(gcdcAcB);
   d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
 
@@ -2913 +3283 @@
 
   if (d == 0)
-  {
-    if (substitute > 1)
-      return N(reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
-  }
+    return N(gcdcAcB);
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
@@ -3018 +3383 @@
 
     if (d0 == 0)
-    {
-      if (substitute > 1)
-        return N(reverseSubst (gcdcAcB, substitute));
-      else
-        return N(gcdcAcB);
-    }
+      return N(gcdcAcB);
     if (d0 >  d)
     {
@@ -3060 +3420 @@
         //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);
-          }
+        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
+        {
+          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);
           return N(gcdcAcB*ppH);
         }
       }
-      else if (fdivides (ppH, A) && fdivides (ppH, B))
-      {
-        if (substitute > 1)
-        {
-          ppH= reverseSubst (ppH, substitute);
-          gcdcAcB= reverseSubst (gcdcAcB, substitute);
-        }
+      else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
         return N(gcdcAcB*ppH);
-      }
     }
     G_m= H;
@@ -3204 +3552 @@
 
         if (d0 == 0)
-        {
-          if (substitute > 1)
-            return N(reverseSubst (gcdcAcB, substitute));
-          else
-            return N(gcdcAcB);
-        }
+          return N(gcdcAcB);
         if (d0 >  d)
         {
@@ -3250 +3593 @@
             //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 (fdivides (ppH, ppA) && fdivides (ppH, ppB))
             {
-              if (substitute > 1)
-              {
-                ppH= reverseSubst (ppH, substitute);
-                gcdcAcB= reverseSubst (gcdcAcB, substitute);
-              }
+              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);
               return N(gcdcAcB*ppH);
             }
           }
-          else if (fdivides (ppH, A) && fdivides (ppH, B))
+          else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
           {
-            if (substitute > 1)
-            {
-              ppH= reverseSubst (ppH, substitute);
-              gcdcAcB= reverseSubst (gcdcAcB, substitute);
-            }
             return N(gcdcAcB*ppH);
           }
@@ -3313 +3646 @@
     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
@@ -3359 +3687 @@
 
   if (d == 0)
-  {
-    if (substitute > 1)
-      return N(reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
-  }
+    return N(gcdcAcB);
+
   d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
 
@@ -3371 +3695 @@
 
   if (d == 0)
-  {
-    if (substitute > 1)
-      return N(reverseSubst (gcdcAcB, substitute));
-    else
-      return N(gcdcAcB);
-  }
+    return N(gcdcAcB);
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
@@ -3519 +3838 @@
 
     if (d0 == 0)
-    {
-      if (substitute > 1)
-        return N(reverseSubst (gcdcAcB, substitute));
-      else
-        return N(gcdcAcB);
-    }
+      return N(gcdcAcB);
     if (d0 >  d)
     {
@@ -3561 +3875 @@
       DEBOUTLN (cerr, "ppH= " << ppH);
 
-      if (fdivides (ppH, A) && fdivides (ppH, B))
-      {
-        if (substitute > 1)
-        {
-          ppH= reverseSubst (ppH, substitute);
-          gcdcAcB= reverseSubst (gcdcAcB, substitute);
-        }
+      if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
         return N(gcdcAcB*ppH);
-      }
     }
     G_m= H;
@@ -3743 +4050 @@
 
         if (d0 == 0)
-        {
-          if (substitute > 1)
-            return N(reverseSubst (gcdcAcB, substitute));
-          else
-            return N(gcdcAcB);
-        }
+          return N(gcdcAcB);
         if (d0 >  d)
         {
@@ -3786 +4088 @@
           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);
-            }
+          if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
             return N(gcdcAcB*ppH);
-          }
         }
 
@@ -4355 +4650 @@
   }
 
-  CanonicalForm cand;
+  CanonicalForm cand, contcand;
   CanonicalForm result;
   int o, t;
@@ -4598 +4893 @@
       if (gcdfound == 1)
       {
-        cand = DD[2] / content( DD[2], Variable(1) );
-        gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
+        contcand= content (DD[2], Variable (1));
+        cand = DD[2] / contcand;
+        if (B_is_F)
+          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
+        else
+          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
 
         if (passToGF && gcdfound)

factory/cf_gcd_smallp.h

@@ -44 +44 @@
                            bool& top_level, CFList& l);
 
+CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G, CanonicalForm& coF, CanonicalForm& coG,
+                           bool& topLevel, CFList& l);
+
 ///GCD of A and B over \f$ F_{p} \f$
 static inline CanonicalForm GCD_small_p (const CanonicalForm& A, const CanonicalForm& B)
@@ -50 +53 @@
   bool top_level= true;
   return GCD_small_p (A, B, top_level, list);
+}
+
+static inline CanonicalForm GCD_small_p (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& coA, CanonicalForm& coB)
+{
+  CFList list;
+  bool top_level= true;
+  return GCD_small_p (A, B, coA, coB, top_level, list);
 }
 
@@ -102 +112 @@
 CFArray
 getMonoms (const CanonicalForm& F);
+
+bool
+terminationTest (const CanonicalForm& F, const CanonicalForm& G,
+                 const CanonicalForm& coF, const CanonicalForm& coG,
+                 const CanonicalForm& cand);
 #endif

factory/fac_ezgcd.cc

@@ -1 @@
-/* emacs edit mode for this file is -*- C++ -*- */
-/* $Id$ */
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR
+\*****************************************************************************/
+/** @file fac_ezgcd.cc
+ *
+ * This file implements the GCD of two multivariate polynomials over Q using
+ * EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor,
+ * Labahnn
+ *
+ * @author Martin Lee
+ *
+ **/
+/*****************************************************************************/
 
 #include "config.h"
@@ -9 +20 @@
 #include "cf_defs.h"
 #include "canonicalform.h"
-#include "fac_util.h"
+#include "fac_util.h" // HEADER for this file
 #include "cf_algorithm.h"
 #include "cf_reval.h"
@@ -16 +27 @@
 #include "fac_distrib.h"
 #include "templates/ftmpl_functions.h"
-
-static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG );
-
-static CanonicalForm ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal );
-
-static CanonicalForm ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound );
-
-static modpk findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG );
-
-static modpk enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk );
+#include "cf_map.h"
+#include "facHensel.h"
+
+static
+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;
+
+  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)
+  {
+    delete [] degsf;
+    delete [] degsg;
+    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
+Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
+                    CFMap & N, const Evaluation& A)
+{
+  int n= F.level();
+  int * degsf= new int [n + 1];
+
+  for (int i = 0; i <= n; i++)
+    degsf[i]= 0;
+
+  degsf= degrees (F, degsf);
+
+  Evaluation result= Evaluation (A.min(), A.max());
+  ASSERT (A.min() == 2, "expected A.min() == 2");
+  ASSERT (A.max() == n, "expected A.max() == n");
+  int max_deg;
+  int k= n;
+  int l= 1;
+  int i= 2;
+  int pos= 2;
+  while (k > 1)
+  {
+    max_deg= degsf [i];
+    while (max_deg == 0)
+    {
+      i++;
+      max_deg= degsf [i];
+    }
+    l= i;
+    for (int j= i + 1; j <=  n; j++)
+    {
+      if (degsf[j] > max_deg)
+      {
+        max_deg= degsf[j];
+        l= j;
+      }
+    }
+
+    if (l <= n)
+    {
+      if (l != pos)
+      {
+        result.setValue (pos, A [l]);
+        M.newpair (Variable (l), Variable (pos));
+        N.newpair (Variable (pos), Variable (l));
+        degsf[l]= 0;
+        l= 2;
+        if (k == 2 && n == 3)
+        {
+          result.setValue (l, A [pos]);
+          M.newpair (Variable (pos), Variable (l));
+          N.newpair (Variable (l), Variable (pos));
+          degsf[pos]= 0;
+        }
+      }
+      else
+      {
+        result.setValue (l, A [l]);
+        degsf [l]= 0;
+      }
+    }
+    pos++;
+    k--;
+    l= 2;
+  }
+
+  delete [] degsf;
+
+  return result;
+}
+
+static inline
+int Hensel (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
+            const CFArray& LeadCoeffs )
+{
+  CFList factors;
+  factors.append (G[1]);
+  factors.append (G[2]);
+
+  CFMap NN, MM;
+  Evaluation A= optimize4Lift (UU, MM, NN, AA);
+
+  CanonicalForm U= MM (UU);
+  CFArray LCs= CFArray (1,2);
+  LCs [1]= MM (LeadCoeffs [1]);
+  LCs [2]= MM (LeadCoeffs [2]);
+
+  CFList evaluation;
+  for (int i= A.min(); i <= A.max(); i++)
+    evaluation.append (A [i]);
+  CFList UEval;
+  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
+
+  if (size (shiftedU)/getNumVars (U) > 500)
+    return -1;
+
+  CFArray shiftedLCs= CFArray (2);
+  CFList shiftedLCsEval1, shiftedLCsEval2;
+  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
+  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
+  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();
+  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
+                        lcs, false);
+
+  for (CFListIterator i= factors; i.hasItem(); i++)
+  {
+    if (!fdivides (i.getItem(), UEval.getFirst()))
+      return 0;
+  }
+
+  int * liftBounds;
+  bool noOneToOne= false;
+  if (U.level() > 2)
+  {
+    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
+    liftBounds[0]= liftBound;
+    for (int i= 1; i < U.level() - 1; i++)
+      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
+    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
+                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
+                                  diophant, noOneToOne);
+    delete [] liftBounds;
+    if (noOneToOne)
+      return 0;
+  }
+  G[1]= factors.getFirst();
+  G[2]= factors.getLast();
+  G[1]= myReverseShift (G[1], evaluation);
+  G[2]= myReverseShift (G[2], evaluation);
+  G[1]= NN (G[1]);
+  G[2]= NN (G[2]);
+  return 1;
+}
+
+static
+bool findeval (const CanonicalForm & F, const CanonicalForm & G,
+               CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
+               REvaluation & b, int delta, int degF, int degG, int maxeval,
+               int & count, int& k, int bound, int& l)
+{
+  if( count == 0 && delta != 0)
+  {
+    if( count++ > maxeval )
+      return false;
+  }
+  if (count > 0)
+  {
+    b.nextpoint(k);
+    if (k == 0)
+      k++;
+    l++;
+    if (l > bound)
+    {
+      l= 1;
+      k++;
+      if (k > tmax (F.level(), G.level()) - 1)
+        return false;
+      b.nextpoint (k);
+    }
+    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
+        {
+          k++;
+          return true;
+        }
+      }
+    }
+    if (k == 0)
+      k++;
+    b.nextpoint(k);
+    l++;
+    if (l > bound)
+    {
+      l= 1;
+      k++;
+      if (k > tmax (F.level(), G.level()) - 1)
+        return false;
+      b.nextpoint (k);
+    }
+    if( count++ > maxeval )
+      return false;
+  }
+}
+
+static CanonicalForm
+ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b,
+        bool internal )
+{
+  bool isRat= isOn (SW_RATIONAL);
+  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
+                lcD, cand, contcand, result;
+  CFArray DD( 1, 2 ), lcDD( 1, 2 );
+  int degF, degG, delta, t, count, maxeval;
+  REvaluation bt;
+  bool gcdfound = false;
+  Variable x = Variable(1);
+  count= 0;
+  maxeval= 200;
+  int o, l;
+  o= 0;
+  l= 1;
+
+  if (!isRat)
+    On (SW_RATIONAL);
+  F= FF*bCommonDen (FF);
+  G= GG*bCommonDen (GG);
+  if (!isRat)
+    Off (SW_RATIONAL);
+
+  CFMap M,N;
+  int smallestDegLev;
+  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
+
+  if (best_level == 0)
+  {
+    DEBDECLEVEL( cerr, "ezgcd" );
+    return G.genOne();
+  }
+
+  F= M (F);
+  G= M (G);
+
+  DEBINCLEVEL( cerr, "ezgcd" );
+  DEBOUTLN( cerr, "FF = " << FF );
+  DEBOUTLN( cerr, "GG = " << GG );
+  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
+  d /= icontent (d);
+  DEBOUTLN( cerr, "f = " << f );
+  DEBOUTLN( cerr, "g = " << g );
+  F /= f; G /= g;
+  if ( F.isUnivariate() && G.isUnivariate() )
+  {
+    DEBDECLEVEL( cerr, "ezgcd" );
+    if(F.mvar()==G.mvar())
+      d*=gcd(F,G);
+    return N (d);
+  }
+
+  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
+  int sizeF= size (F);
+  int sizeG= size (G);
+
+
+  if (!isRat)
+    On (SW_RATIONAL);
+  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
+  {
+    Off(SW_USE_EZGCD);
+    result=gcd( F, G );
+    On(SW_USE_EZGCD);
+    if (!isRat)
+      Off (SW_RATIONAL);
+    result /= icontent (result);
+    DEBDECLEVEL( cerr, "ezgcd" );
+    return N (d*result);
+  }
+
+  if ( gcd_test_one( F, G, false ) )
+  {
+    DEBDECLEVEL( cerr, "ezgcd" );
+    if (!isRat)
+      Off (SW_RATIONAL);
+    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 );
+  t = tmax( F.level(), G.level() );
+  if ( ! internal )
+    b = REvaluation( 2, t, IntRandom( 25 ) );
+  while ( ! gcdfound )
+  {
+    /// ---> A2
+    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
+    DEBOUTLN( cerr, "F = " << F );
+    DEBOUTLN( cerr, "G = " << G );
+    if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
+                   o, 25, l))
+    {
+      Off(SW_USE_EZGCD);
+      result=gcd( F, G );
+      On(SW_USE_EZGCD);
+      if (!isRat)
+        Off (SW_RATIONAL);
+      DEBDECLEVEL( cerr, "ezgcd" );
+      result /= icontent (result);
+      return N (d*result);
+    }
+    DEBOUTLN( cerr, "found evaluation b = " << b );
+    DEBOUTLN( cerr, "F(b) = " << Fb );
+    DEBOUTLN( cerr, "G(b) = " << Gb );
+    DEBOUTLN( cerr, "D(b) = " << Db );
+    delta = degree( Db );
+    /// ---> A3
+    if ( delta == 0 )
+    {
+      DEBDECLEVEL( cerr, "ezgcd" );
+      if (!isRat)
+        Off (SW_RATIONAL);
+      return N (d);
+    }
+    /// ---> A4
+    //deltaold = delta;
+    while ( 1 ) 
+    {
+      bt = b;
+      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
+                     o, 25,l ))
+      {
+        Off(SW_USE_EZGCD);
+        result=gcd( F, G );
+        On(SW_USE_EZGCD);
+        if (!isRat)
+          Off (SW_RATIONAL);
+        DEBDECLEVEL( cerr, "ezgcd" );
+        result /= icontent (result);
+        return N (d*result);
+      }
+      int dd=degree( Dbt );
+      if ( dd /*degree( Dbt )*/ == 0 )
+      {
+        DEBDECLEVEL( cerr, "ezgcd" );
+        if (!isRat)
+          Off (SW_RATIONAL);
+        return N (d);
+      }
+      if ( dd /*degree( Dbt )*/ == delta )
+        break;
+      else  if ( dd /*degree( Dbt )*/ < delta )
+      {
+        delta = dd /*degree( Dbt )*/;
+        b = bt;
+        Db = Dbt; Fb = Fbt; Gb = Gbt;
+      }
+    }
+    DEBOUTLN( cerr, "now after A4, delta = " << delta );
+    /// ---> A5
+    if (delta == degF)
+    {
+      if (degF <= degG  && fdivides (F, G))
+      {
+        DEBDECLEVEL( cerr, "ezgcd" );
+        if (!isRat)
+          Off (SW_RATIONAL);
+        return N (d*F);
+      }
+      else
+        delta--;
+    }
+    else if (delta == degG)
+    {
+      if (degG <= degF && fdivides( G, F ))
+      {
+        DEBDECLEVEL( cerr, "ezgcd" );
+        if (!isRat)
+          Off (SW_RATIONAL);
+        return N (d*G);
+      }
+      else
+        delta--;
+    }
+    if ( delta != degF && delta != degG )
+    {
+      /// ---> A6
+      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
+        Off(SW_USE_EZGCD);
+        result=gcd( F, G );
+        On(SW_USE_EZGCD);
+        if (!isRat)
+          Off (SW_RATIONAL);
+        DEBDECLEVEL( cerr, "ezgcd" );
+        result /= icontent (result);
+        return N (d*result);
+      }
+      /// ---> A7
+      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
+      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
+      DEBOUTLN( cerr, "(hensel) B    = " << B );
+      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
+      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
+      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
+      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
+      gcdfound= Hensel (B*lcD, DD, b, lcDD);
+      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
+
+      if (gcdfound == -1)
+      {
+        Off (SW_USE_EZGCD);
+        result= gcd (F,G);
+        On (SW_USE_EZGCD);
+        if (!isRat)
+          Off (SW_RATIONAL);
+        DEBDECLEVEL( cerr, "ezgcd" );
+        result /= icontent (result);
+        return N (d*result);
+      }
+
+      if (gcdfound)
+      {
+        contcand= content (DD[2], Variable (1));
+        cand = DD[2] / contcand;
+        if (B_is_F)
+          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
+        else
+          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
+      }
+      /// ---> A8 (gcdfound)
+    }
+    delta--;
+  }
+  /// ---> A9
+  DEBDECLEVEL( cerr, "ezgcd" );
+  cand *= bCommonDen (cand);
+  if (!isRat)
+    Off (SW_RATIONAL);
+  cand /= icontent (cand);
+  return N (d*cand);
+}
 
 CanonicalForm
 ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
 {
-    REvaluation b;
-    return ezgcd( FF, GG, b, false );
-}
-
-static CanonicalForm
-ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal )
-{
-    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, t;
-    REvaluation bt;
-    bool gcdfound = false;
-    Variable x = Variable(1);
-    modpk bound;
-
-    DEBINCLEVEL( cerr, "ezgcd" );
-    DEBOUTLN( cerr, "FF = " << FF );
-    DEBOUTLN( cerr, "GG = " << GG );
-    f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
-    DEBOUTLN( cerr, "f = " << f );
-    DEBOUTLN( cerr, "g = " << g );
-    F = FF / f; G = GG / g;
-    if ( F.isUnivariate() && G.isUnivariate() )
-    {
-        DEBDECLEVEL( cerr, "ezgcd" );
-        if(F.mvar()==G.mvar())
-          d*=gcd(F,G);
-        return d;
-    }
-    else  if ( gcd_test_one( F, G, false ) )
-    {
-        DEBDECLEVEL( cerr, "ezgcd" );
-        return d;
-    }
-    lcF = LC( F, x ); lcG = LC( G, x );
-    lcD = gcd( lcF, lcG );
-    delta = 0;
-    degF = degree( F, x ); degG = degree( G, x );
-    t = tmax( F.level(), G.level() );
-    bound = findBound( F, G, lcF, lcG, degF, degG );
-    if ( ! internal )
-        b = REvaluation( 2, t, IntRandom( 25 ) );
-    while ( ! gcdfound ) {
-        /// ---> A2
-        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
-        DEBOUTLN( cerr, "F = " << F );
-        DEBOUTLN( cerr, "G = " << G );
-        findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
-        DEBOUTLN( cerr, "found evaluation b = " << b );
-        DEBOUTLN( cerr, "F(b) = " << Fb );
-        DEBOUTLN( cerr, "G(b) = " << Gb );
-        DEBOUTLN( cerr, "D(b) = " << Db );
-        delta = degree( Db );
-        /// ---> A3
-        if ( delta == 0 )
-        {
-          DEBDECLEVEL( cerr, "ezgcd" );
-          return d;
-        }
-        /// ---> A4
-        //deltaold = delta;
-        while ( 1 ) {
-            bt = b;
-            findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
-            int dd=degree( Dbt );
-            if ( dd /*degree( Dbt )*/ == 0 )
-            {
-                DEBDECLEVEL( cerr, "ezgcd" );
-                return d;
-            }
-            if ( dd /*degree( Dbt )*/ == delta )
-                break;
-            else  if ( dd /*degree( Dbt )*/ < delta ) {
-                delta = dd /*degree( Dbt )*/;
-                b = bt;
-                Db = Dbt; Fb = Fbt; Gb = Gbt;
-            }
-        }
-        DEBOUTLN( cerr, "now after A4, delta = " << delta );
-        /// ---> A5
-        if ( degF <= degG && delta == degF && fdivides( F, G ) )
-        {
-            DEBDECLEVEL( cerr, "ezgcd" );
-            return d*F;
-        }
-        if ( degG < degF && delta == degG && fdivides( G, F ) )
-        {
-            DEBDECLEVEL( cerr, "ezgcd" );
-            return d*G;
-        }
-        if ( delta != degF && delta != degG ) {
-            bool B_is_F;
-            /// ---> A6
-            CanonicalForm xxx;
-            //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
-            DD[1] = Fb / Db;
-            xxx= gcd( DD[1], Db );
-            DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
-            DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
-            DEBOUTLN( cerr, "Db = " << Db );
-            if (xxx.inCoeffDomain()) {
-                B = F;
-                DD[2] = Db;
-                lcDD[1] = lcF;
-                lcDD[2] = lcF;
-                B *= lcF;
-                B_is_F=true;
-            }
-            //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
-            else
-            {
-              DD[1] = Gb / Db;
-              xxx=gcd( DD[1], Db );
-              DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
-              DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
-              DEBOUTLN( cerr, "Db = " << Db );
-              if (xxx.inCoeffDomain())
-              {
-                B = G;
-                DD[2] = Db;
-                lcDD[1] = lcG;
-                lcDD[2] = lcG;
-                B *= lcG;
-                B_is_F=false;
-              }
-              else
-              {
-#ifdef DEBUGOUTPUT
-                CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
-                DEBDECLEVEL( cerr, "ezgcd" );
-                return dummyres;
-#else
-                return d * ezgcd_specialcase( F, G, b, bound );
-#endif
-              }
-            }
-            /// ---> A7
-            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
-            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
-            bound = enlargeBound( B, DD[2], DD[1], bound );
-            DEBOUTLN( cerr, "(hensel) B    = " << B );
-            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
-            DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
-            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
-            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
-            gcdfound = Hensel( B, DD, lcDD, b, bound, x );
-            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
-
-            if (gcdfound)
-            {
-              CanonicalForm xxx=content(DD[2],Variable(1));
-              CanonicalForm cand=DD[2] / xxx; //content(DD[2],Variable(1));
-#if 0
-              gcdfound= fdivides(cand,G) &&  fdivides(cand,F);
-#else
-              if (B_is_F /*B==F*lcF*/)
-              {
-                DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" -> " << G%cand );
-                gcdfound= fdivides(cand,G);
-              }
-              else
-              {
-                DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" -> " << F%cand);
-                gcdfound= fdivides(cand,F);
-              }
-#endif
-            }
-            /// ---> A8 (gcdfound)
-        }
-        delta++;
-    }
-    /// ---> A9
-    DEBDECLEVEL( cerr, "ezgcd" );
-    return d * DD[2] / content(DD[2],Variable(1));
-}
-
-static CanonicalForm
-ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & /*b*/, const modpk & /*bound*/ )
-{
-    CanonicalForm d;
-#if 1
-    Off(SW_USE_EZGCD);
-    //bool ntl0=isOn(SW_USE_NTL_GCD_0);
-    //Off(SW_USE_NTL_GCD_0);
-    //bool ntlp=isOn(SW_USE_NTL_GCD_P);
-    //Off(SW_USE_NTL_GCD_P);
-    d=gcd( F, G );
-    //if (ntl0) On(SW_USE_NTL_GCD_0);
-    //if (ntlp) On(SW_USE_NTL_GCD_P);
-    On(SW_USE_EZGCD);
-    return d;
-#else
-    /*
-     * I am not sure, if these lines are much of use.
-     * Remove?
-     */
-    DEBOUTLN( cerr, "ezgcd: special case" );
-    CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt;
-    CFArray DD( 1, 2 ), lcDD( 1, 2 );
-    Variable x = Variable( 1 );
-    bool gcdfound;
-
-    Dt = 1;
-    /// ---> S1
-    DEBOUTLN( cerr, "ezgcdspec: (S1)" );
-    Ft = ezgcd( F, F.deriv( x ) );
-    if ( Ft.isOne() )
-    {
-        // In this case F is squarefree and we came here by bad chance
-        // (means: bad evaluation point).  Try again with another
-        // evaluation point.  Bug fix (?) by JS.  The bad example was
-        // gcd.debug -ocr /+USE_EZGCD/@12/CB
-        //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)'
-        //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
-        b.nextpoint();
-        return ezgcd( F, G, b, true );
-    }
-
-    DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft );
-    L = F / Ft;
-    /// ---> S2
-    DEBOUTLN( cerr, "ezgcdspec: (S2)" );
-
-    L = ezgcd( L, G, b, true );
-    DEBOUTLN( cerr, "ezgcdspec: (S2) done, Ds = " << Ds );
-    Gt = G / L;
-    Ds = L; Dt = L;
-    Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
-    d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
-    while ( true ) {
-        /// ---> S3
-        DEBOUTLN( cerr, "ezgcdspec: (S3)" );
-        DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
-        DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
-        DD[1] = gcd( Lb, gcd( Fb, Gb ) );
-        /// ---> S4
-        DEBOUTLN( cerr, "ezgcdspec: (S4)" );
-        if ( degree( DD[1] ) == 0 )
-            return Ds;
-        /// ---> S5
-        DEBOUTLN( cerr, "ezgcdspec: (S5)" );
-        DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
-        DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
-        Db = DD[1];
-        if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
-            DEBOUTLN( cerr, "ezgcdspec: (S55)" );
-            lcDD[2] = LC( L, x );
-            lcDD[1] = d;
-            DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
-            DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
-            LL = L * d;
-            modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
-            DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
-            DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
-            DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
-            DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
-            DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
-            DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
-            DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
-            DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
-            gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
-            ASSERT( gcdfound, "fatal error in algorithm" );
-            Dt = pp( DD[1] );
-        }
-        DEBOUTLN( cerr, "ezgcdspec: (S7)" );
-        Ds = Ds * Dt;
-        Fb = Fb / Db;
-        Gb = Gb / Db;
-    }
-    // this point is never reached, but to avoid compiler warnings let's return a value
-    return 0;
-#endif
-}
-
-static void
-findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG )
-{
-    bool ok;
-    if ( delta != 0 )
-        b.nextpoint();
-    DEBOUTLN( cerr, "ezgcd: (findeval) F = " << F  <<", G="<< G);
-    DEBOUTLN( cerr, "ezgcd: (findeval) degF = " << degF << ", degG="<<degG );
-    do {
-        DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
-        Fb = b( F );
-        ok = degree( Fb ) == degF;
-        if ( ok ) {
-            Gb = b( G );
-            ok = degree( Gb ) == degG;
-        }
-
-        if ( ok )
-        {
-            Db = gcd( Fb, Gb );
-            if ( delta > 0 )
-              ok = degree( Db ) < delta;
-        }
-        if ( ! ok )
-        {
-            b.nextpoint();
-        }
-    } while ( ! ok );
-}
-
-static modpk
-enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk )
-{
-    DEBOUTLN( cerr, "ezgcd: (enlarge bound) F      = " << F );
-    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb     = " << Lb );
-    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db     = " << Db );
-    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb % p = " << mod( Lb, pk.getp() ) );
-    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db % p = " << mod( Db, pk.getp() ) );
-
-    CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) *
-        tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) );
-    int p = pk.getp();
-    int k = pk.getk();
-    CanonicalForm bound = pk.getpk();
-    while ( bound < limit ) {
-        k++;
-        bound *= p;
-    }
-    k *= 2;
-    DEBOUTLN( cerr, "ezgcd: (enlarge bound) newbound = " << power( CanonicalForm( p ), k ) );
-    return modpk( p, k );
-}
-
-static modpk
-findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG )
-{
-    CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) *
-        gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) );
-    int p, i = 0;
-    do {
-        p = cf_getBigPrime( i );
-        i++;
-    } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() );
-    CanonicalForm bound = p;
-    i = 1;
-    while ( bound < limit ) {
-        i++;
-        bound *= p;
-    }
-    return modpk( p, i );
-}
+  REvaluation b;
+  return ezgcd( FF, GG, b, false );
+}
+