Changeset dec1024 in git

Timestamp:

Sep 20, 2010, 5:17:36 PM (13 years ago)

Author:

Martin Lee <martinlee84@…>

Branches:

(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')

Children:

a0b5e27a5fb475f020a02a82976ffe87d3c097ae

Parents:

3f46965e1ad849ccb6d39d9e57cca2ecbbf10d86

Message:

added heuristics to modular gcd


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

File:

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

factory/cf_gcd_smallp.cc

@@ -50 +50 @@
 static inline
 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
-                CFMap & N, bool& top_level) 
+                CFMap & N, bool& topLevel) 
 { 
   int n= tmax (F.level(), G.level());
@@ -67 +67 @@
   int both_zero= 0;
 
-  if (top_level) 
+  if (topLevel) 
   {
     for (int i= 1; i <= n; i++) 
@@ -88 +88 @@
     }
 
-    if (both_non_zero == 0) return 0;
+    if (both_non_zero == 0) 
+    {
+      delete [] degsf;
+      delete [] degsg;
+      return 0;
+    }
 
     // map Variables which do not occur in both polynomials to higher levels
@@ -201 +206 @@
 }
 
+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 seems to be 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());   
+}
+
+static inline CanonicalForm 
+uni_content (const CanonicalForm & F);
+
+CanonicalForm
+uni_content (const CanonicalForm& F, const Variable& x)
+{
+  if (F.inCoeffDomain())
+    return F.genOne();
+  if (F.level() == x.level() && F.isUnivariate())
+    return F;
+  if (F.level() != x.level() && F.isUnivariate())
+    return F.genOne();
+    
+  if (x.level() != 1)
+  {
+    CanonicalForm f= swapvar (F, x, Variable (1));
+    CanonicalForm result= uni_content (f);
+    return swapvar (result, x, Variable (1));
+  }
+  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$ 
@@ -231 +376 @@
     return c;
   }
+}
+
+CanonicalForm 
+extractContents (const CanonicalForm& F, const CanonicalForm& G, 
+                 CanonicalForm& contentF, CanonicalForm& contentG, 
+                 CanonicalForm& ppF, CanonicalForm& ppG, const int d)
+{
+  CanonicalForm uniContentF, uniContentG, gcdcFcG;
+  contentF= 1;
+  contentG= 1;
+  ppF= F;
+  ppG= G;
+  CanonicalForm result= 1;
+  for (int i= 1; i <= d; i++)
+  {
+    uniContentF= uni_content (F, Variable (i));
+    uniContentG= uni_content (G, Variable (i));
+    gcdcFcG= gcd (uniContentF, uniContentG);
+    contentF *= uniContentF;
+    contentG *= uniContentG;
+    ppF /= uniContentF;
+    ppG /= uniContentG;
+    result *= gcdcFcG;
+  }
+  return result;
 }
 
@@ -336 +506 @@
 
 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 
-/// l and top_level are only used internally, output is monic
+/// 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& top_level) 
+                  Variable & alpha, CFList& l, bool& topLevel) 
 { 
   CanonicalForm A= F;
@@ -354 +524 @@
  
   CFMap M,N;
-  int best_level= myCompress (A, B, M, N, top_level);
+  int best_level= myCompress (A, B, M, N, topLevel);
 
   if (best_level == 0) return B.genOne();
@@ -367 +537 @@
     return N (gcd(A,B)); 
   
+  int substitute= substituteCheck (A, B);
+  
+  if (substitute > 1)
+    subst (A, B, A, B, substitute);
+
   CanonicalForm cA, cB;    // content of A and B
   CanonicalForm ppA, ppB;    // primitive part of A and B
   CanonicalForm gcdcAcB;
 
-  cA = uni_content (A);
-  cB = uni_content (B);
-
-  gcdcAcB= gcd (cA, cB);
-
-  ppA= A/cA;
-  ppB= B/cB;
+  if (topLevel)
+  {
+    if (best_level <= 2)
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
+    else 
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);  
+  }
+  else
+  {
+    cA = uni_content (A);
+    cB = uni_content (B); 
+    gcdcAcB= gcd (cA, cB);
+    ppA= A/cA;
+    ppB= B/cB;
+  }
 
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
@@ -385 +568 @@
   lcB= uni_lcoeff (ppB);
   
-  if (fdivides (lcA, lcB)) { 
+  if (fdivides (lcA, lcB)) 
+  { 
     if (fdivides (A, B))
       return F/Lc(F);    
@@ -399 +583 @@
   int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
 
-  if (d == 0) return N(gcdcAcB);
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N (reverseSubst (gcdcAcB, substitute));
+    else 
+      return N(gcdcAcB);
+  }
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
   if (d0 < d) 
     d= d0; 
-  if (d == 0) return N(gcdcAcB);
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N (reverseSubst (gcdcAcB, substitute));
+    else 
+      return N(gcdcAcB);
+  }
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
@@ -413 +609 @@
   H= 0;
   bool fail= false;
-  top_level= false;
+  topLevel= false;
   bool inextension= false;
   Variable V_buf= alpha;
@@ -461 +657 @@
       G_random_element= 
       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
-                        list, top_level);
+                        list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -472 +668 @@
       G_random_element= 
       GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
-                        list, top_level);
+                        list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -480 +676 @@
     d0= totaldegree (G_random_element, Variable(2), 
                      Variable (G_random_element.level()));
-    if (d0 == 0) return N(gcdcAcB); 
+    if (d0 == 0)
+    {
+      if (substitute > 1)
+        return N (reverseSubst (gcdcAcB, substitute));
+      else 
+        return N(gcdcAcB);
+    }
     if (d0 >  d) 
     {
@@ -522 +724 @@
         DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
         if (fdivides (ppH, A) && fdivides (ppH, B))         
+        {
+          if (substitute > 1)
+          {
+            ppH= reverseSubst (ppH, substitute);
+            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+          }
           return N(gcdcAcB*ppH);
+        }
       }
       else if (fdivides (ppH, A) && fdivides (ppH, B)) 
-        return N(gcdcAcB*(ppH/Lc(ppH)));
+      {
+        if (substitute > 1)
+        {
+          ppH= reverseSubst (ppH, substitute);
+          gcdcAcB= reverseSubst (gcdcAcB, substitute);
+        }
+        return N(gcdcAcB*ppH);
+      }
     }
 
@@ -580 +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& top_level) 
+        CFList& l, bool& topLevel) 
 { 
   CanonicalForm A= F;
@@ -593 +809 @@
  
   CFMap M,N;
-  int best_level= myCompress (A, B, M, N, top_level);
+  int best_level= myCompress (A, B, M, N, topLevel);
 
   if (best_level == 0) return B.genOne();
@@ -606 +822 @@
     return N (gcd (A, B)); 
 
+  int substitute= substituteCheck (A, B);
+  
+  if (substitute > 1)
+    subst (A, B, A, B, substitute);
+
   CanonicalForm cA, cB;    // content of A and B
   CanonicalForm ppA, ppB;    // primitive part of A and B
   CanonicalForm gcdcAcB;
 
-  cA = uni_content (A);
-  cB = uni_content (B);
-
-  gcdcAcB= gcd (cA, cB);
-
-  ppA= A/cA;
-  ppB= B/cB;
+  if (topLevel)
+  {
+    if (best_level <= 2)
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
+    else 
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);  
+  }
+  else
+  {
+    cA = uni_content (A);
+    cB = uni_content (B); 
+    gcdcAcB= gcd (cA, cB);
+    ppA= A/cA;
+    ppB= B/cB;
+  }
 
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
@@ -638 +867 @@
 
   int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
-  if (d == 0) return N(gcdcAcB);
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N (reverseSubst (gcdcAcB, substitute));
+    else 
+      return N(gcdcAcB);
+  }
   int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
   if (d0 < d) 
     d= d0; 
-  if (d == 0) return N(gcdcAcB);
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N (reverseSubst (gcdcAcB, substitute));
+    else 
+      return N(gcdcAcB);
+  }
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
@@ -652 +893 @@
   H= 0;
   bool fail= false;
-  top_level= false;
+  topLevel= false;
   bool inextension= false;
   int p;
@@ -694 +935 @@
       TIMING_START (gcd_recursion);
       G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
-                                list, top_level);
+                                list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -704 +945 @@
       TIMING_START (gcd_recursion);
       G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
-                                list, top_level);
+                                list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -717 +958 @@
       {
         setCharacteristic (p, k, gf_name_buf);
-        return N(gcdcAcB); 
+        {
+          if (substitute > 1)
+            return N (reverseSubst (gcdcAcB, substitute));
+          else
+            return N(gcdcAcB); 
+        } 
       }
       else
-        return N(gcdcAcB);
+      {
+        if (substitute > 1)
+          return N (reverseSubst (gcdcAcB, substitute));
+        else
+          return N(gcdcAcB); 
+      }
     } 
     if (d0 >  d) 
@@ -759 +1010 @@
           ppH= GFMapDown (ppH, k);
           DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
+          if (substitute > 1)
+          {
+            ppH= reverseSubst (ppH, substitute);
+            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+          }
           setCharacteristic (p, k, gf_name_buf);
           return N(gcdcAcB*ppH);
@@ -766 +1022 @@
       {
         if (fdivides (ppH, A) && fdivides (ppH, B)) 
+        {
+          if (substitute > 1)
+          {
+            ppH= reverseSubst (ppH, substitute);
+            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+          }
           return N(gcdcAcB*ppH);
+        }
       }
     }
@@ -829 +1092 @@
 
 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G, 
-                           bool& top_level, CFList& l) 
+                           bool& topLevel, CFList& l) 
 {
   CanonicalForm A= F;
@@ -842 +1105 @@
 
   CFMap M,N;
-  int best_level= myCompress (A, B, M, N, top_level);
+  int best_level= myCompress (A, B, M, N, topLevel);
 
   if (best_level == 0) return B.genOne();
@@ -855 +1118 @@
     return N (gcd (A, B));
 
+  int substitute= substituteCheck (A, B);
+  
+  if (substitute > 1)
+    subst (A, B, A, B, substitute);
+
   CanonicalForm cA, cB;    // content of A and B
   CanonicalForm ppA, ppB;    // primitive part of A and B
   CanonicalForm gcdcAcB;
-  cA = uni_content (A);
-  cB = uni_content (B); 
-  gcdcAcB= gcd (cA, cB);
-  ppA= A/cA;
-  ppB= B/cB;
+
+  if (topLevel)
+  {
+    if (best_level <= 2)
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
+    else 
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);  
+  }
+  else
+  {
+    cA = uni_content (A);
+    cB = uni_content (B); 
+    gcdcAcB= gcd (cA, cB);
+    ppA= A/cA;
+    ppB= B/cB;
+  }
 
   CanonicalForm lcA, lcB;  // leading coefficients of A and B
@@ -885 +1164 @@
   int d0;
 
-  if (d == 0) return N(gcdcAcB);
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N (reverseSubst (gcdcAcB, substitute));
+    else 
+      return N(gcdcAcB);
+  }
   d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
 
@@ -891 +1176 @@
     d= d0;
 
-  if (d == 0) return N(gcdcAcB);
+  if (d == 0) 
+  {
+    if (substitute > 1)
+      return N (reverseSubst (gcdcAcB, substitute));
+    else 
+      return N(gcdcAcB);
+  }
 
   CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
@@ -899 +1190 @@
   G_m= 0;
   Variable alpha, V_buf;
-  int p, k;
   bool fail= false;
   bool inextension= false;
   bool inextensionextension= false;
-  top_level= false;
-  CanonicalForm CF_buf;
+  topLevel= false;
   CFList source, dest;
-  CanonicalForm gcdcheck;
   do 
   {
@@ -919 +1207 @@
       TIMING_START (gcd_recursion);
       G_random_element= 
-      GCD_small_p (ppA (random_element,x), ppB (random_element,x), top_level, 
+      GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 
       list);
       TIMING_END_AND_PRINT (gcd_recursion, 
@@ -931 +1219 @@
       G_random_element= 
       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
-                        list, top_level);
+                        list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -942 +1230 @@
       CFList list;
       CanonicalForm mipo;
-      int deg= 3;
+      int deg= 2;
       do {
         mipo= randomIrredpoly (deg, x); 
@@ -955 +1243 @@
       G_random_element= 
       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
-                        list, top_level);
+                        list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -1001 +1289 @@
       G_random_element= 
       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
-                        list, top_level);
+                        list, topLevel);
       TIMING_END_AND_PRINT (gcd_recursion, 
                             "time for recursive call: ");
@@ -1013 +1301 @@
     if (d0 == 0)
     {
-      return N(gcdcAcB); 
+      if (substitute > 1)
+        return N (reverseSubst (gcdcAcB, substitute));
+      else
+        return N(gcdcAcB); 
     }
     if (d0 >  d) 
@@ -1050 +1341 @@
       DEBOUTLN (cerr, "ppH= " << ppH);
       if (fdivides (ppH, A) && fdivides (ppH, B)) 
+      {
+        if (substitute > 1)
+        {
+          ppH= reverseSubst (ppH, substitute);
+          gcdcAcB= reverseSubst (gcdcAcB, substitute);
+        }
         return N(gcdcAcB*ppH);
+      }
     }