Changeset 806c18 in git for factory/cf_gcd_smallp.cc

Timestamp:

Nov 15, 2010, 4:34:57 PM (13 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')

Children:

7c3bca08c96331a56864c1d35b8c2e8ff2e0be89

Parents:

c840d97af622b4e4da8761738b540e21144f716b

Message:

format

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

File:

• factory/cf_gcd_smallp.cc (modified) (87 diffs)

factory/cf_gcd_smallp.cc

@@ -6 +6 @@
  * @date 22.10.2009
  *
- * This file implements the GCD of two polynomials over \f$ F_{p} \f$ , 
- * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms 
+ * This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
+ * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms
  * for Computer Algebra" by Geddes, Czapor, Labahnn
  *
@@ -13 +13 @@
  *   (c) by The SINGULAR Team, see LICENSE file
  *
- * @internal 
+ * @internal
  * @version \$Id$
  *
@@ -46 +46 @@
 TIMING_DEFINE_PRINT(newton_interpolation);
 
-/// compressing two polynomials F and G, M is used for compressing, 
+/// compressing two polynomials F and G, M is used for compressing,
 /// N to reverse the compression
 static inline
 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
-                CFMap & N, bool& topLevel) 
-{ 
+                CFMap & N, bool& topLevel)
+{
   int n= tmax (F.level(), G.level());
   int * degsf= new int [n + 1];
@@ -58 +58 @@
   for (int i = 0; i <= n; i++)
     degsf[i]= degsg[i]= 0;
-   
+
   degsf= degrees (F, degsf);
   degsg= degrees (G, degsg);
-  
+
   int both_non_zero= 0;
   int f_zero= 0;
@@ -67 +67 @@
   int both_zero= 0;
 
-  if (topLevel) 
-  {
-    for (int i= 1; i <= n; i++) 
-    {
-      if (degsf[i] != 0 && degsg[i] != 0) 
+  if (topLevel)
+  {
+    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()) 
+      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
       {
         f_zero++;
         continue;
       }
-      if (degsg[i] == 0 && degsf[i] && i <= F.level()) 
+      if (degsg[i] == 0 && degsf[i] && i <= F.level())
       {
         g_zero++;
@@ -88 +88 @@
     }
 
-    if (both_non_zero == 0) 
+    if (both_non_zero == 0)
     {
       delete [] degsf;
@@ -98 +98 @@
     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) 
+    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));
@@ -109 +109 @@
         k++;
       }
-      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 
-      {
-        if (l + g_zero + both_non_zero != i) 
+      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));
@@ -119 +119 @@
       }
     }
- 
+
     // sort Variables x_{i} in increasing order of
-    // min(deg_{x_{i}}(f),deg_{x_{i}}(g)) 
+    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
     int m= tmin (F.level(), G.level());
     int max_min_deg;
@@ -127 +127 @@
     l= 0;
     int i= 1;
-    while (k > 0) 
+    while (k > 0)
     {
       max_min_deg= tmin (degsf[i], degsg[i]);
-      while (max_min_deg == 0) 
+      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) 
+      for (int j= i + 1; j <=  m; j++)
+      {
+        if (tmin (degsf[j],degsg[j]) >= max_min_deg)
         {
           max_min_deg= tmin (degsf[j], degsg[j]);
-          l= j; 
+          l= j;
         }
       }
-      if (l != 0) 
-      {
-        if (l != k) 
+      if (l != 0)
+      {
+        if (l != k)
         {
           M.newpair (Variable (l), Variable(k));
@@ -153 +153 @@
           l= 0;
         }
-        else 
+        else
         {
           degsf[l]= 0;
@@ -159 +159 @@
           l= 0;
         }
-      }  
-      else if (l == 0) 
-      {
-        if (i != k) 
+      }
+      else if (l == 0)
+      {
+        if (i != k)
         {
           M.newpair (Variable (i), Variable (k));
@@ -169 +169 @@
           degsg[i]= 0;
         }
-        else 
+        else
         {
           degsf[i]= 0;
@@ -175 +175 @@
         }
         i++;
-      } 
+      }
       k--;
     }
   }
-  else 
+  else
   {
     //arrange Variables such that no gaps occur
-    for (int i= 1; i <= n; i++) 
-    {
-      if (degsf[i] == 0 && degsg[i] == 0) 
+    for (int i= 1; i <= n; i++)
+    {
+      if (degsf[i] == 0 && degsg[i] == 0)
       {
         both_zero++;
         continue;
       }
-      else 
-      {
-        if (both_zero != 0) 
+      else
+      {
+        if (both_zero != 0)
         {
           M.newpair (Variable (i), Variable (i - both_zero));
@@ -207 +207 @@
 
 int
-substituteCheck (const CanonicalForm& F, const CanonicalForm& G) 
+substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
 {
   if (F.inCoeffDomain() || G.inCoeffDomain())
@@ -214 +214 @@
   if (degree (F, x) <= 1 || degree (G, x) <= 1)
     return 0;
-  CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive 
-  CanonicalForm g= swapvar (G, G.mvar(), x); 
+  CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive
+  CanonicalForm g= swapvar (G, G.mvar(), x);
   int sizef= 0;
-  int sizeg= 0; 
+  int sizeg= 0;
   for (CFIterator i= f; i.hasTerms(); i++, sizef++)
   {
@@ -240 +240 @@
     expg [j]= i.exp();
   }
-  
+
   int indf= sizef - 1;
   int indg= sizeg - 1;
@@ -247 +247 @@
   if (expg[indg] == 0)
     indg--;
-    
+
   if (expg[indg] != expf [indf] || (expg[indg] == 1 && expf[indf] == 1))
   {
@@ -265 +265 @@
     }
   }
-  
+
   for (int i= indg - 1; i >= 0; i--)
   {
@@ -282 +282 @@
 
 // substiution
-void 
-subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A, 
+void
+subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A,
        CanonicalForm& B, const int d
       )
@@ -293 +293 @@
     return;
   }
-  
+
   CanonicalForm C= 0;
-  CanonicalForm D= 0;  
+  CanonicalForm D= 0;
   Variable x= Variable (1);
   CanonicalForm f= swapvar (F, x, F.mvar());
@@ -307 +307 @@
 }
 
-CanonicalForm 
+CanonicalForm
 reverseSubst (const CanonicalForm& F, const int d)
 {
   if (d == 1)
     return F;
-  
+
   Variable x= Variable (1);
   if (degree (F, x) <= 0)
@@ -320 +320 @@
   for (CFIterator i= f; i.hasTerms(); i++)
     result += i.coeff()*power (f.mvar(), d*i.exp());
-  return swapvar (result, x, F.mvar());   
-}
-
-static inline CanonicalForm 
+  return swapvar (result, x, F.mvar());
+}
+
+static inline CanonicalForm
 uni_content (const CanonicalForm & F);
 
@@ -335 +335 @@
   if (F.level() != x.level() && F.isUnivariate())
     return F.genOne();
-    
+
   if (x.level() != 1)
   {
@@ -344 +344 @@
   else
     return uni_content (F);
-} 
-
-/// compute the content of F, where F is considered as an element of 
-/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$ 
-static inline CanonicalForm 
-uni_content (const CanonicalForm & F) 
-{  
+}
+
+/// compute the content of F, where F is considered as an element of
+/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
+static inline CanonicalForm
+uni_content (const CanonicalForm & F)
+{
   if (F.inBaseDomain())
     return F.genOne();
@@ -360 +360 @@
 
   int l= F.level();
-  if (l == 2) 
+  if (l == 2)
     return content(F);
-  else 
+  else
   {
     CanonicalForm pol, c = 0;
     CFIterator i = F;
-    for (; i.hasTerms(); i++) 
-    {
-      pol= i.coeff(); 
+    for (; i.hasTerms(); i++)
+    {
+      pol= i.coeff();
       pol= uni_content (pol);
       c= gcd (c, pol);
@@ -378 +378 @@
 }
 
-CanonicalForm 
-extractContents (const CanonicalForm& F, const CanonicalForm& G, 
-                 CanonicalForm& contentF, CanonicalForm& contentG, 
+CanonicalForm
+extractContents (const CanonicalForm& F, const CanonicalForm& G,
+                 CanonicalForm& contentF, CanonicalForm& contentG,
                  CanonicalForm& ppF, CanonicalForm& ppG, const int d)
 {
@@ -407 +407 @@
 /// is dp.
 static inline
-CanonicalForm uni_lcoeff (const CanonicalForm& F) 
+CanonicalForm uni_lcoeff (const CanonicalForm& F)
 {
   if (F.level() <= 1)
-    return F; 
+    return F;
   else
   {
@@ -418 +418 @@
     {
       if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
-        return uni_lcoeff (i.coeff()); 
+        return uni_lcoeff (i.coeff());
     }
   }
@@ -441 +441 @@
 }
 
-/// compute a random element a of \f$ F_{p}(\alpha ) \f$ , 
+/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
 /// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
-/// fail if there are no field elements left which have not been used before 
-static inline CanonicalForm 
+/// fail if there are no field elements left which have not been used before
+static inline CanonicalForm
 randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
-               bool & fail) 
+               bool & fail)
 {
   fail= false;
@@ -457 +457 @@
   int d= degree (mipo);
   int bound= ipower (p, d);
-  do 
+  do
   {
     if (list.length() == bound)
@@ -464 +464 @@
       break;
     }
-    if (list.length() < p) 
+    if (list.length() < p)
     {
       random= genFF.generate();
@@ -470 +470 @@
         random= genFF.generate();
     }
-    else 
+    else
     {
       random= genAlgExt.generate();
@@ -476 +476 @@
         random= genAlgExt.generate();
     }
-    if (F (random, x) == 0) 
+    if (F (random, x) == 0)
     {
       list.append (random);
@@ -485 +485 @@
 }
 
-/// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$ 
+/// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$
 /// we do not extend \f$ F_{p}(\alpha ) \f$ itself but extend \f$ F_{p} \f$ ,
-/// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$ 
+/// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$
 static inline
-void choose_extension (const int& d, const int& num_vars, Variable& beta) 
+void choose_extension (const int& d, const int& num_vars, Variable& beta)
 {
   int p= getCharacteristic();
@@ -496 +496 @@
   ZZ_pX NTLirredpoly;
   //TODO: replace d by max_{i} (deg_x{i}(f))
-  int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p)); 
+  int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p));
   int m= degree (getMipo (beta));
   if (i <= 1)
     i= 2;
-  BuildIrred (NTLirredpoly, i*m); 
-  CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1)); 
-  beta= rootOf (mipo); 
-}
-
-/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 
+  BuildIrred (NTLirredpoly, i*m);
+  CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1));
+  beta= rootOf (mipo);
+}
+
+/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
 /// l and topLevel are only used internally, output is monic
 /// based on Alg. 7.2. as described in "Algorithms for
 /// Computer Algebra" by Geddes, Czapor, Labahn
-CanonicalForm 
-GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 
-                  Variable & alpha, CFList& l, bool& topLevel) 
-{ 
+CanonicalForm
+GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
+                  Variable & alpha, CFList& l, bool& topLevel)
+{
   CanonicalForm A= F;
   CanonicalForm B= G;
@@ -522 +522 @@
   if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
   if (F == G) return F/Lc(F);
- 
+
   CFMap M,N;
   int best_level= myCompress (A, B, M, N, topLevel);
@@ -533 +533 @@
   Variable x= Variable(1);
 
-  //univariate case 
-  if (A.isUnivariate() && B.isUnivariate()) 
-    return N (gcd(A,B)); 
-  
+  //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);
@@ -550 +550 @@
     if (best_level <= 2)
       gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
-    else 
-      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);  
+    else
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
   }
   else
   {
     cA = uni_content (A);
-    cB = uni_content (B); 
+    cB = uni_content (B);
     gcdcAcB= gcd (cA, cB);
     ppA= A/cA;
@@ -563 +563 @@
 
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
-  CanonicalForm gcdlcAlcB; 
+  CanonicalForm gcdlcAlcB;
 
   lcA= uni_lcoeff (ppA);
   lcB= uni_lcoeff (ppB);
-  
-  if (fdivides (lcA, lcB)) 
-  { 
+
+  if (fdivides (lcA, lcB))
+  {
     if (fdivides (A, B))
-      return F/Lc(F);    
+      return F/Lc(F);
   }
   if (fdivides (lcB, lcA))
-  { 
-    if (fdivides (B, A)) 
+  {
+    if (fdivides (B, A))
       return G/Lc(G);
   }
@@ -587 +587 @@
     if (substitute > 1)
       return N (reverseSubst (gcdcAcB, substitute));
-    else 
+    else
       return N(gcdcAcB);
   }
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
-  if (d0 < d) 
-    d= d0; 
+  if (d0 < d)
+    d= d0;
   if (d == 0)
   {
     if (substitute > 1)
       return N (reverseSubst (gcdcAcB, substitute));
-    else 
+    else
       return N(gcdcAcB);
   }
@@ -614 +614 @@
   CanonicalForm prim_elem, im_prim_elem;
   CFList source, dest;
-  do 
+  do
   {
     random_element= randomElement (m, V_buf, l, fail);
-    if (fail) 
+    if (fail)
     {
       source= CFList();
@@ -638 +638 @@
       DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
       inextension= true;
-      for (CFListIterator i= l; i.hasItem(); i++) 
-        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 
+      for (CFListIterator i= l; i.hasItem(); i++)
+        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
                              im_prim_elem, source, dest);
       m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
-      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 
+      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
                           source, dest);
       ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
-      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 
+      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
                          source, dest);
 
@@ -655 +655 @@
       CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= 
+      G_random_element=
       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
                         list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     }
-    else 
+    else
     {
       CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= 
+      G_random_element=
       GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
                         list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     }
 
-    d0= totaldegree (G_random_element, Variable(2), 
+    d0= totaldegree (G_random_element, Variable(2),
                      Variable (G_random_element.level()));
     if (d0 == 0)
@@ -680 +680 @@
       if (substitute > 1)
         return N (reverseSubst (gcdcAcB, substitute));
-      else 
+      else
         return N(gcdcAcB);
     }
-    if (d0 >  d) 
+    if (d0 >  d)
     {
       if (!find (l, random_element))
@@ -690 +690 @@
     }
 
-    G_random_element= 
+    G_random_element=
     (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
     * G_random_element;
 
-    d0= totaldegree (G_random_element, Variable(2), 
+    d0= totaldegree (G_random_element, Variable(2),
                      Variable(G_random_element.level()));
-    if (d0 <  d) 
+    if (d0 <  d)
     {
       m= gcdlcAlcB;
@@ -706 +706 @@
     TIMING_START (newton_interpolation);
     H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
-    TIMING_END_AND_PRINT (newton_interpolation, 
+    TIMING_END_AND_PRINT (newton_interpolation,
                           "time for newton interpolation: ");
 
-    //termination test 
-    if (uni_lcoeff (H) == gcdlcAlcB) 
+    //termination test
+    if (uni_lcoeff (H) == gcdlcAlcB)
     {
       cH= uni_content (H);
       ppH= H/cH;
-      if (inextension) 
-      {
-        CFList u, v; 
+      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= 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, A) && fdivides (ppH, B))
         {
           if (substitute > 1)
@@ -733 +733 @@
         }
       }
-      else if (fdivides (ppH, A) && fdivides (ppH, B)) 
+      else if (fdivides (ppH, A) && fdivides (ppH, B))
       {
         if (substitute > 1)
@@ -752 +752 @@
 }
 
-/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a 
+/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
 /// univariate polynomial, returns fail if there are no field elements left
 /// which have not been used before
@@ -766 +766 @@
   int d= getGFDegree();
   int bound= ipower (p, d);
-  do 
+  do
   {
     if (list.length() == bound)
@@ -775 +775 @@
     if (list.length() < 1)
       random= 0;
-    else 
+    else
     {
       random= genGF.generate();
@@ -781 +781 @@
         random= genGF.generate();
     }
-    if (F (random, x) == 0) 
+    if (F (random, x) == 0)
     {
       list.append (random);
@@ -796 +796 @@
 /// the size of the base field is bounded by 2^16, output is monic
 CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
-        CFList& l, bool& topLevel) 
-{ 
+        CFList& l, bool& topLevel)
+{
   CanonicalForm A= F;
   CanonicalForm B= G;
@@ -807 +807 @@
   if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
   if (F == G) return F/Lc(F);
- 
+
   CFMap M,N;
   int best_level= myCompress (A, B, M, N, topLevel);
@@ -818 +818 @@
   Variable x= Variable(1);
 
-  //univariate case 
-  if (A.isUnivariate() && B.isUnivariate()) 
-    return N (gcd (A, B)); 
+  //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);
@@ -835 +835 @@
     if (best_level <= 2)
       gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
-    else 
-      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);  
+    else
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
   }
   else
   {
     cA = uni_content (A);
-    cB = uni_content (B); 
+    cB = uni_content (B);
     gcdcAcB= gcd (cA, cB);
     ppA= A/cA;
@@ -848 +848 @@
 
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
-  CanonicalForm gcdlcAlcB; 
+  CanonicalForm gcdlcAlcB;
 
   lcA= uni_lcoeff (ppA);
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB)) 
-  { 
+  if (fdivides (lcA, lcB))
+  {
     if (fdivides (A, B))
-      return F; 
-  }   
-  if (fdivides (lcB, lcA)) 
-  { 
-    if (fdivides (B, A)) 
+      return F;
+  }
+  if (fdivides (lcB, lcA))
+  {
+    if (fdivides (B, A))
       return G;
   }
@@ -871 +871 @@
     if (substitute > 1)
       return N (reverseSubst (gcdcAcB, substitute));
-    else 
+    else
       return N(gcdcAcB);
   }
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
-  if (d0 < d) 
-    d= d0; 
+  if (d0 < d)
+    d= d0;
   if (d == 0)
   {
     if (substitute > 1)
       return N (reverseSubst (gcdcAcB, substitute));
-    else 
+    else
       return N(gcdcAcB);
   }
@@ -900 +900 @@
   int expon;
   char gf_name_buf= gf_name;
-  do 
+  do
   {
     random_element= GFRandomElement (m, l, fail);
-    if (fail) 
-    { 
+    if (fail)
+    {
       int num_vars= tmin (getNumVars(A), getNumVars(B));
       num_vars--;
@@ -911 +911 @@
       if (expon < 2)
         expon= 2;
-      kk= getGFDegree(); 
-      if (ipower (p, kk*expon) < (1 << 16)) 
+      kk= getGFDegree();
+      if (ipower (p, kk*expon) < (1 << 16))
         setCharacteristic (p, kk*(int)expon, 'b');
-      else 
+      else
       {
         expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
@@ -922 +922 @@
       inextension= true;
       fail= false;
-      for (CFListIterator i= l; i.hasItem(); i++) 
+      for (CFListIterator i= l; i.hasItem(); i++)
         i.getItem()= GFMapUp (i.getItem(), kk);
       m= GFMapUp (m, kk);
@@ -934 +934 @@
       CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
+      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
                                 list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     }
-    else 
+    else
     {
       CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
+      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
                                 list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     }
 
-    d0= totaldegree (G_random_element, Variable(2), 
+    d0= totaldegree (G_random_element, Variable(2),
                      Variable (G_random_element.level()));
-    if (d0 == 0) 
-    {
-      if (inextension) 
+    if (d0 == 0)
+    {
+      if (inextension)
       {
         setCharacteristic (p, k, gf_name_buf);
@@ -962 +962 @@
             return N (reverseSubst (gcdcAcB, substitute));
           else
-            return N(gcdcAcB); 
-        } 
+            return N(gcdcAcB);
+        }
       }
       else
@@ -970 +970 @@
           return N (reverseSubst (gcdcAcB, substitute));
         else
-          return N(gcdcAcB); 
-      }
-    } 
-    if (d0 >  d) 
+          return N(gcdcAcB);
+      }
+    }
+    if (d0 >  d)
     {
       if (!find (l, random_element))
@@ -980 +980 @@
     }
 
-    G_random_element= 
+    G_random_element=
     (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
      G_random_element;
-    d0= totaldegree (G_random_element, Variable(2), 
+    d0= totaldegree (G_random_element, Variable(2),
                      Variable (G_random_element.level()));
 
-    if (d0 < d) 
+    if (d0 < d)
     {
       m= gcdlcAlcB;
@@ -998 +998 @@
     TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: ");
 
-    //termination test 
-    if (uni_lcoeff (H) == gcdlcAlcB) 
+    //termination test
+    if (uni_lcoeff (H) == gcdlcAlcB)
     {
       cH= uni_content (H);
       ppH= H/cH;
-      if (inextension) 
-      {
-        if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k))) 
+      if (inextension)
+      {
+        if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k)))
         {
           DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
@@ -1019 +1019 @@
         }
       }
-      else 
-      {
-        if (fdivides (ppH, A) && fdivides (ppH, B)) 
+      else
+      {
+        if (fdivides (ppH, A) && fdivides (ppH, B))
         {
           if (substitute > 1)
@@ -1042 +1042 @@
 
 /// F is assumed to be an univariate polynomial in x,
-/// computes a random monic irreducible univariate polynomial of random 
+/// computes a random monic irreducible univariate polynomial of random
 /// degree < i in x which does not divide F
-CanonicalForm 
-randomIrredpoly (int i, const Variable & x) 
+CanonicalForm
+randomIrredpoly (int i, const Variable & x)
 {
   int p= getCharacteristic();
   ZZ NTLp= to_ZZ (p);
   ZZ_p::init (NTLp);
-  ZZ_pX NTLirredpoly; 
+  ZZ_pX NTLirredpoly;
   CanonicalForm CFirredpoly;
   BuildIrred (NTLirredpoly, i + 1);
   CFirredpoly= convertNTLZZpX2CF (NTLirredpoly, x);
   return CFirredpoly;
-} 
+}
 
 static inline
@@ -1067 +1067 @@
   int p= getCharacteristic();
   int bound= p;
-  do 
+  do
   {
     if (list.length() == bound)
@@ -1076 +1076 @@
     if (list.length() < 1)
       random= 0;
-    else 
+    else
     {
       random= genFF.generate();
@@ -1082 +1082 @@
         random= genFF.generate();
     }
-    if (F (random, x) == 0) 
+    if (F (random, x) == 0)
     {
       list.append (random);
@@ -1091 +1091 @@
 }
 
-CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G, 
-                           bool& topLevel, CFList& l) 
+CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
+                           bool& topLevel, CFList& l)
 {
   CanonicalForm A= F;
@@ -1114 +1114 @@
   Variable x= Variable (1);
 
-  //univariate case 
-  if (A.isUnivariate() && B.isUnivariate()) 
+  //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);
@@ -1131 +1131 @@
     if (best_level <= 2)
       gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
-    else 
-      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);  
+    else
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
   }
   else
   {
     cA = uni_content (A);
-    cB = uni_content (B); 
+    cB = uni_content (B);
     gcdcAcB= gcd (cA, cB);
     ppA= A/cA;
@@ -1144 +1144 @@
 
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
-  CanonicalForm gcdlcAlcB; 
+  CanonicalForm gcdlcAlcB;
   lcA= uni_lcoeff (ppA);
   lcB= uni_lcoeff (ppB);
 
-  if (fdivides (lcA, lcB)) 
-  { 
+  if (fdivides (lcA, lcB))
+  {
     if (fdivides (A, B))
       return F/Lc(F);
-  }    
+  }
   if (fdivides (lcB, lcA))
-  { 
-    if (fdivides (B, A)) 
+  {
+    if (fdivides (B, A))
       return G/Lc(G);
   }
 
-  gcdlcAlcB= gcd (lcA, lcB); 
-  
+  gcdlcAlcB= gcd (lcA, lcB);
+
   int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
   int d0;
@@ -1168 +1168 @@
     if (substitute > 1)
       return N (reverseSubst (gcdcAcB, substitute));
-    else 
+    else
       return N(gcdcAcB);
   }
   d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
 
-  if (d0 < d) 
+  if (d0 < d)
     d= d0;
 
-  if (d == 0) 
+  if (d == 0)
   {
     if (substitute > 1)
       return N (reverseSubst (gcdcAcB, substitute));
-    else 
+    else
       return N(gcdcAcB);
   }
@@ -1195 +1195 @@
   topLevel= false;
   CFList source, dest;
-  do 
+  do
   {
     if (inextension)
@@ -1206 +1206 @@
       CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= 
-      GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 
+      G_random_element=
+      GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel,
       list);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
@@ -1215 +1215 @@
     else if (!fail && inextension)
     {
-      CFList list; 
+      CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= 
-      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
+      G_random_element=
+      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
                         list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
-      DEBOUTLN (cerr, "G_random_element= " << G_random_element);    
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     }
     else if (fail && !inextension)
@@ -1232 +1232 @@
       int deg= 2;
       do {
-        mipo= randomIrredpoly (deg, x); 
+        mipo= randomIrredpoly (deg, x);
         alpha= rootOf (mipo);
         inextension= true;
         fail= false;
-        random_element= randomElement (m, alpha, l, fail); 
+        random_element= randomElement (m, alpha, l, fail);
         deg++;
-      } while (fail); 
+      } while (fail);
       list= CFList();
       TIMING_START (gcd_recursion);
-      G_random_element= 
-      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
+      G_random_element=
+      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
                         list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
@@ -1260 +1260 @@
       CanonicalForm prim_elem, im_prim_elem;
       prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
-       
+
       ASSERT (!prim_fail, "failure in integer factorizer");
       if (prim_fail)
@@ -1271 +1271 @@
 
       inextensionextension= true;
-      for (CFListIterator i= l; i.hasItem(); i++) 
-        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 
+      for (CFListIterator i= l; i.hasItem(); i++)
+        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
                              im_prim_elem, source, dest);
       m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
-      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 
+      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
                           source, dest);
       ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
       ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
-      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 
+      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
                          source, dest);
       fail= false;
@@ -1287 +1287 @@
       CFList list;
       TIMING_START (gcd_recursion);
-      G_random_element= 
+      G_random_element=
       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
                         list, topLevel);
-      TIMING_END_AND_PRINT (gcd_recursion, 
+      TIMING_END_AND_PRINT (gcd_recursion,
                             "time for recursive call: ");
       DEBOUTLN (cerr, "G_random_element= " << G_random_element);
       alpha= V_buf;
-    } 
-
-    d0= totaldegree (G_random_element, Variable(2), 
+    }
+
+    d0= totaldegree (G_random_element, Variable(2),
                      Variable (G_random_element.level()));
 
@@ -1304 +1304 @@
         return N (reverseSubst (gcdcAcB, substitute));
       else
-        return N(gcdcAcB); 
-    }
-    if (d0 >  d) 
-    { 
+        return N(gcdcAcB);
+    }
+    if (d0 >  d)
+    {
       if (!find (l, random_element))
         l.append (random_element);
@@ -1314 +1314 @@
 
     G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
-                       *G_random_element;  
-
-    
-    d0= totaldegree (G_random_element, Variable(2), 
+                       *G_random_element;
+
+
+    d0= totaldegree (G_random_element, Variable(2),
                      Variable(G_random_element.level()));
 
-    if (d0 <  d) 
+    if (d0 <  d)
     {
       m= gcdlcAlcB;
@@ -1327 +1327 @@
       d= d0;
     }
-     
+
     TIMING_START (newton_interpolation);
     H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
-    TIMING_END_AND_PRINT (newton_interpolation, 
+    TIMING_END_AND_PRINT (newton_interpolation,
                           "time for newton_interpolation: ");
 
     //termination test
-    if (uni_lcoeff (H) == gcdlcAlcB) 
+    if (uni_lcoeff (H) == gcdlcAlcB)
     {
       cH= uni_content (H);
@@ -1340 +1340 @@
       ppH /= Lc (ppH);
       DEBOUTLN (cerr, "ppH= " << ppH);
-      if (fdivides (ppH, A) && fdivides (ppH, B)) 
+      if (fdivides (ppH, A) && fdivides (ppH, B))
       {
         if (substitute > 1)