Changeset fe2d4c in git for factory/algext.cc

Timestamp:

Jun 20, 2011, 2:10:50 PM (13 years ago)

Author:

Martin Lee <martinlee84@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

4e35a894dd1ecea25a42d6b12a182338771e5732

Parents:

35eb6cf319007843f3388d3f071f3bee96188fbd

Message:

added better compression of polynomials to tryBrownGCD
added better termination test to QGCD


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

File:

• factory/algext.cc (modified) (8 diffs)

factory/algext.cc

@@ -14 +14 @@
 #endif
 
+#include "templates/ftmpl_functions.h"
 #include "cf_defs.h"
 #include "canonicalform.h"
@@ -23 +24 @@
 #include "cf_map.h"
 #include "cf_generator.h"
+
+/// compressing two polynomials F and G, M is used for compressing,
+/// N to reverse the compression
+static
+int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
+                CFMap & N, bool topLevel)
+{
+  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);
+
+  int both_non_zero= 0;
+  int f_zero= 0;
+  int g_zero= 0;
+  int both_zero= 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())
+      {
+        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 increasing order of
+    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
+    int m= tmax (F.level(), G.level());
+    int min_max_deg;
+    k= both_non_zero;
+    l= 0;
+    int i= 1;
+    while (k > 0)
+    {
+      if (degsf [i] != 0 && degsg [i] != 0)
+        min_max_deg= tmax (degsf[i], degsg[i]);
+      else
+        min_max_deg= 0;
+      while (min_max_deg == 0)
+      {
+        i++;
+        min_max_deg= tmax (degsf[i], degsg[i]);
+        if (degsf [i] != 0 && degsg [i] != 0)
+          min_max_deg= tmax (degsf[i], degsg[i]);
+        else
+          min_max_deg= 0;
+      }
+      for (int j= i + 1; j <=  m; j++)
+      {
+        if (tmax (degsf[j],degsg[j]) <= min_max_deg && degsf[j] != 0 && degsg [j] != 0)
+        {
+          min_max_deg= tmax (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--;
+    }
+  }
+  else
+  {
+    //arrange Variables such that no gaps occur
+    for (int i= 1; i <= n; i++)
+    {
+      if (degsf[i] == 0 && degsg[i] == 0)
+      {
+        both_zero++;
+        continue;
+      }
+      else
+      {
+        if (both_zero != 0)
+        {
+          M.newpair (Variable (i), Variable (i - both_zero));
+          N.newpair (Variable (i - both_zero), Variable (i));
+        }
+      }
+    }
+  }
+
+  delete [] degsf;
+  delete [] degsg;
+
+  return 1;
+}
 
 CanonicalForm reduce(const CanonicalForm & f, const CanonicalForm & M)
@@ -199 +367 @@
 
 
-void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail )
+void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
 { // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
   // M is assumed to be monic
@@ -261 +429 @@
   }
   CFMap MM,NN;
-  CFArray ps(1,2);
-  ps[1] = F;
-  ps[2] = G;
-  compress(ps,MM,NN); // maps MM, NN are created
+  int lev= myCompress (F, G, MM, NN, topLevel);
+  if (lev == 0)
+  {
+    result= 1;
+    return;
+  }
   CanonicalForm f=MM(F);
   CanonicalForm g=MM(G);
@@ -347 +517 @@
     if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
       continue;
-    tryBrownGCD( f(alpha, Variable(1)), g(alpha, Variable(1)), M, g_image, fail ); // recursive call with one var less
+    tryBrownGCD( f(alpha, Variable(1)), g(alpha, Variable(1)), M, g_image, fail, false ); // recursive call with one var less
     if(fail)
       return;
@@ -488 +658 @@
   q = 1;
   D = 0;
+  CanonicalForm test= 0;
+  bool equal= false;
   for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
   {
@@ -530 +702 @@
       setReduce(a,true); // reduce expressions modulo mipo
       On( SW_RATIONAL ); // needed by fdivides
-      if( fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
+      if (test != tmp)
+        test= tmp;
+      else
+        equal= true; // modular image did not add any new information
+      if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
       {
         Off( SW_RATIONAL );
@@ -762 +938 @@
   while(true)
   {
-    tryInvert( Lc(result), M, inv, fail );
-    if(fail)
-      return;
-    // here: Lc(result) is invertible modulo M
-    q = Lc(P)*inv*power( x, P.degree(x)-result.degree(x) );
-    rem = reduce( P - q*result, M );
+    tryDivrem (P, result, q, rem, inv, M, fail);
+    if(fail)
+      return;
     if( rem.isZero() )
     {