Index: factory/cf_gcd.cc
===================================================================
--- factory/cf_gcd.cc	(revision cb7827a842fd2e9878498f940e9f7735fdbe488e)
+++ factory/cf_gcd.cc	(revision 59778389f5075553bef793320610ce298d74b18d)
@@ -1238,7 +1238,7 @@
   CanonicalForm Dn, test= 0;
   cl =  gcd (f.lc(),g.lc());
+  CanonicalForm gcdcfcg= gcd (cf, cg);
   CanonicalForm b= 1;
   int minCommonDeg= 0;
-  CanonicalForm gcdcfcg= gcd (cf, cg);
   for (i= tmax (f.level(), g.level()); i > 0; i--)
   {
@@ -1271,5 +1271,5 @@
     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 );
@@ -1278,9 +1278,6 @@
     //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 = GCD_small_p (mapinto (f), mapinto (g), cofp, cogp);
     Dp /=Dp.lc();
-    cofp= mapinto (f)/Dp;
-    cogp= mapinto (g)/Dp;
     cofp /= lc (cofp);
     cogp /= lc (cogp);
@@ -1348,5 +1345,6 @@
         equal= true;
       //Dn /=vcontent(Dn,Variable(1));
-      if ((equal || q > b) || ((cofn*Dn==f) && (cogn*Dn==g)))
+      if (((abs (cofn*Dn)==abs (f)) && (abs (cogn*Dn)==abs (g))) ||
+          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
       {
         //printf(" -> success\n");
Index: factory/cf_gcd_smallp.cc
===================================================================
--- factory/cf_gcd_smallp.cc	(revision cb7827a842fd2e9878498f940e9f7735fdbe488e)
+++ factory/cf_gcd_smallp.cc	(revision 59778389f5075553bef793320610ce298d74b18d)
@@ -425,4 +425,19 @@
 }
 
+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
@@ -431,20 +446,57 @@
 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);
@@ -455,5 +507,10 @@
   //univariate case
   if (A.isUnivariate() && B.isUnivariate())
-    return N (gcd(A,B));
+  {
+    CanonicalForm result= gcd (A, B);
+    coF= N (A/result);
+    coG= N (B/result);
+    return N (result);
+  }
 
   CanonicalForm cA, cB;    // content of A and B
@@ -484,4 +541,6 @@
   if (compressConvexDense)
   {
+    CanonicalForm bufcA= cA;
+    CanonicalForm bufcB= cB;
     cA= content (ppA, 1);
     cB= content (ppB, 1);
@@ -489,10 +548,21 @@
     ppB /= cB;
     gcdcAcB *= gcd (cA, cB);
+    cA *= bufcA;
+    cB *= bufcB;
     if (ppA.isUnivariate() || ppB.isUnivariate())
     {
       if (ppA.level() == ppB.level())
-        return N (gcd (ppA, ppB)*gcdcAcB);
+      {
+        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);
+      }
     }
 
@@ -511,7 +581,16 @@
     {
       if (ppA.level() == ppB.level())
-        return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
+      {
+        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);
+      }
     }
   }
@@ -523,5 +602,5 @@
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB))
+  /*if (fdivides (lcA, lcB))
   {
     if (fdivides (A, B))
@@ -532,5 +611,5 @@
     if (fdivides (B, A))
       return G/Lc(G);
-  }
+  }*/
 
   gcdlcAlcB= gcd (lcA, lcB);
@@ -539,5 +618,9 @@
 
   if (d == 0)
+  {
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
     return N(gcdcAcB);
+  }
 
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
@@ -545,12 +628,19 @@
     d= d0;
   if (d == 0)
+  {
+    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;
@@ -560,9 +650,9 @@
   CanonicalForm prim_elem, im_prim_elem;
   CFList source, dest;
-  int bound= tmin (degree (ppA, 1), degree (ppB, 1));
-  int count= 0;
+  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)
     {
@@ -579,5 +669,7 @@
       {
         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);
@@ -586,4 +678,6 @@
         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,
@@ -605,4 +699,6 @@
       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);
@@ -611,12 +707,15 @@
       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,
@@ -629,5 +728,6 @@
       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,
@@ -643,5 +743,9 @@
 
     if (d0 == 0)
+    {
+      coF= N (ppA*(cA/gcdcAcB));
+      coG= N (ppB*(cB/gcdcAcB));
       return N(gcdcAcB);
+    }
     if (d0 >  d)
     {
@@ -654,4 +758,8 @@
     (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())
@@ -667,42 +775,68 @@
       G_m= 0;
       d= d0;
-      count= 0;
+      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: ");
 
-    count++;
     //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 ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB)))
-        {
+        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+             (ppA == ppCoF*ppH && ppB == 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= 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 ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB)))
+      else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+                (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+                (ppA == ppCoF*ppH && ppB == 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);
       }
@@ -710,4 +844,6 @@
 
     G_m= H;
+    coF_m= coF;
+    coG_m= coG;
     newtonPoly= newtonPoly*(x - random_element);
     m= m*(x - random_element);
@@ -755,4 +891,18 @@
 }
 
+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
@@ -760,21 +910,59 @@
 /// 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);
@@ -785,5 +973,10 @@
   //univariate case
   if (A.isUnivariate() && B.isUnivariate())
-    return N (gcd (A, B));
+  {
+    CanonicalForm result= gcd (A, B);
+    coF= N (A/result);
+    coG= N (B/result);
+    return N (result);
+  }
 
   CanonicalForm cA, cB;    // content of A and B
@@ -814,4 +1007,6 @@
   if (compressConvexDense)
   {
+    CanonicalForm bufcA= cA;
+    CanonicalForm bufcB= cB;
     cA= content (ppA, 1);
     cB= content (ppB, 1);
@@ -819,10 +1014,21 @@
     ppB /= cB;
     gcdcAcB *= gcd (cA, cB);
+    cA *= bufcA;
+    cB *= bufcB;
     if (ppA.isUnivariate() || ppB.isUnivariate())
     {
       if (ppA.level() == ppB.level())
-        return N (gcd (ppA, ppB)*gcdcAcB);
+      {
+        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);
+      }
     }
 
@@ -841,7 +1047,16 @@
     {
       if (ppA.level() == ppB.level())
-        return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
+      {
+        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);
+      }
     }
   }
@@ -853,14 +1068,26 @@
   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);
@@ -868,20 +1095,31 @@
   int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
   if (d == 0)
+  {
+    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)
     d= d0;
   if (d == 0)
+  {
+    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;
@@ -890,9 +1128,9 @@
   int expon;
   char gf_name_buf= gf_name;
-  int bound= tmin (degree (ppA, 1), degree (ppB, 1));
-  int count= 0;
+  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)
     {
@@ -915,12 +1153,17 @@
       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,
@@ -933,4 +1176,5 @@
       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,
@@ -949,4 +1193,6 @@
       if (inextension)
         setCharacteristic (p, k, gf_name_buf);
+      coF= N (ppA*(cA/gcdcAcB));
+      coG= N (ppB*(cB/gcdcAcB));
       return N(gcdcAcB);
     }
@@ -961,4 +1207,10 @@
     (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),
@@ -973,26 +1225,50 @@
       G_m= 0;
       d= d0;
-      count= 0;
+      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: ");
-
-    count++;
+    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 ((count == bound) || (fdivides(ppH, ppA) && fdivides(ppH, ppB)))
+        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+             (ppA == ppCoF*ppH && ppB == 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 (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);
           setCharacteristic (p, k, gf_name_buf);
           return N(gcdcAcB*ppH);
@@ -1001,8 +1277,17 @@
       else
       {
-        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
+      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+           (ppA == ppCoF*ppH && ppB == 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);
         }
@@ -1011,4 +1296,6 @@
 
     G_m= H;
+    coF_m= coF;
+    coG_m= coG;
     newtonPoly= newtonPoly*(x - random_element);
     m= m*(x - random_element);
@@ -1068,21 +1355,72 @@
 }
 
-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);
@@ -1093,5 +1431,10 @@
   //univariate case
   if (A.isUnivariate() && B.isUnivariate())
-    return N (gcd (A, B));
+  {
+    CanonicalForm result= gcd (A, B);
+    coF= N (A/result);
+    coG= N (B/result);
+    return N (result);
+  }
 
   CanonicalForm cA, cB;    // content of A and B
@@ -1122,4 +1465,6 @@
   if (compressConvexDense)
   {
+    CanonicalForm bufcA= cA;
+    CanonicalForm bufcB= cB;
     cA= content (ppA, 1);
     cB= content (ppB, 1);
@@ -1127,10 +1472,21 @@
     ppB /= cB;
     gcdcAcB *= gcd (cA, cB);
+    cA *= bufcA;
+    cB *= bufcB;
     if (ppA.isUnivariate() || ppB.isUnivariate())
     {
       if (ppA.level() == ppB.level())
-        return N (gcd (ppA, ppB)*gcdcAcB);
+      {
+        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);
+      }
     }
 
@@ -1149,7 +1505,16 @@
     {
       if (ppA.level() == ppB.level())
-        return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
+      {
+        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);
+      }
     }
   }
@@ -1160,5 +1525,5 @@
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB))
+  /*if (fdivides (lcA, lcB))
   {
     if (fdivides (A, B))
@@ -1169,5 +1534,5 @@
     if (fdivides (B, A))
       return G/Lc(G);
-  }
+  }*/
 
   gcdlcAlcB= gcd (lcA, lcB);
@@ -1177,5 +1542,9 @@
 
   if (d == 0)
+  {
+    coF= N (ppA*(cA/gcdcAcB));
+    coG= N (ppB*(cB/gcdcAcB));
     return N(gcdcAcB);
+  }
 
   d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
@@ -1185,10 +1554,19 @@
 
   if (d == 0)
+  {
+    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;
@@ -1196,14 +1574,14 @@
   bool inextension= false;
   bool inextensionextension= false;
-  topLevel= false;
+  bool topLevel2= false;
   CFList source, dest;
-  int bound= tmin (degree (ppA, 1), degree (ppB, 1));
-  int count= 0;
+  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)
@@ -1212,6 +1590,7 @@
       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: ");
@@ -1223,6 +1602,7 @@
       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: ");
@@ -1241,5 +1621,5 @@
         inextension= true;
         fail= false;
-        random_element= randomElement (m, alpha, l, fail);
+        random_element= randomElement (m*lcA*lcB, alpha, l, fail);
         deg++;
       } while (fail);
@@ -1248,6 +1628,7 @@
       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: ");
@@ -1269,5 +1650,7 @@
       {
         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);
@@ -1276,4 +1659,6 @@
         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,
@@ -1296,4 +1681,6 @@
       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);
@@ -1302,12 +1689,15 @@
       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: ");
@@ -1322,5 +1712,10 @@
 
     if (d0 == 0)
+    {
+      coF= N (ppA*(cA/gcdcAcB));
+      coG= N (ppB*(cB/gcdcAcB));
       return N(gcdcAcB);
+    }
+
     if (d0 >  d)
     {
@@ -1333,4 +1728,8 @@
                        *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())
@@ -1346,24 +1745,43 @@
       G_m= 0;
       d= d0;
-      count= 0;
+      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: ");
 
-    count++;
     //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 ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB)))
+      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
+           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
+           (ppA == ppCoF*ppH && ppB == 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);
       }
@@ -1371,4 +1789,6 @@
 
     G_m= H;
+    coF_m= coF;
+    coG_m= coG;
     newtonPoly= newtonPoly*(x - random_element);
     m= m*(x - random_element);
@@ -4452,5 +4872,5 @@
       {
         cand = DD[2] / content( DD[2], Variable(1) );
-        gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
+        gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); //TODO use multiplication instead of fdivides
 
         if (passToGF && gcdfound)
Index: factory/cf_gcd_smallp.h
===================================================================
--- factory/cf_gcd_smallp.h	(revision cb7827a842fd2e9878498f940e9f7735fdbe488e)
+++ factory/cf_gcd_smallp.h	(revision 59778389f5075553bef793320610ce298d74b18d)
@@ -44,4 +44,7 @@
                            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,4 +53,11 @@
   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);
 }
 
