Changeset f613a24 in git

Timestamp:

Jun 13, 2012, 9:31:39 AM (12 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

02d8afb153e54b4dca7a1c7e782b73b7df40b090

Parents:

ae1e243d742c878b2d7dcb278a6785540ef2b894

Message:

fix: Singular/minpoly.{cc|h}

bugfixes and optimizations
on behalf of Sebastian Jambor
version sent in via e-mail on June 6, 2012, 17:00:31 MESZ for master

Location:

Singular

Files:

• minpoly.cc (modified) (17 diffs)
• minpoly.h (modified) (4 diffs)

Singular/minpoly.cc

@@ -1 @@
-/***********************************************
- * Copyright (C) 2011 Sebastian Jambor         *
- * sebastian@momo.math.rwth-aachen.de          *
- ***********************************************/
+/*************************************************
+ * Author: Sebastian Jambor, 2011                *
+ * GPL (e-mail from June 6, 2012, 17:00:31 MESZ) *
+ * sebastian@momo.math.rwth-aachen.de            *
+ ************************************************/
 
 
@@ -8 +9 @@
 #include "config.h"
 #include<kernel/mod2.h>
+
 //#include<iomanip>
 
@@ -67 +69 @@
       for(int j = piv; j < n + rows + 1; j++)
       {
-        unsigned long tmp = multMod (matrix[i][j], x, p);
-        tmp = p - tmp;
-        tmprow[j] += tmp;
-        tmprow[j] %= p;
+        if (matrix[i][j] != 0)
+        {
+          unsigned long tmp = multMod (matrix[i][j], x, p);
+          tmp = p - tmp;
+          tmprow[j] += tmp;
+          if (tmprow[j] >= p)
+          {
+            tmprow[j] -= p;
+          }
+        }
       }
     }
@@ -182 +190 @@
 
   pivots = new unsigned[n];
+
+  nonPivots = new unsigned[n];
+
+  for (int i = 0; i < n; i++)
+  {
+     nonPivots[i] = i;
+  }
+
   rows = 0;
 }
@@ -187 +203 @@
 NewVectorMatrix::~NewVectorMatrix ()
 {
+  delete nonPivots;
   delete pivots;
 
@@ -226 +243 @@
     {
       // subtract row[i] times the i-th row
-      for(int j = piv; j < n; j++)
-      {
-        unsigned long tmp = multMod (matrix[i][j], x, p);
-        tmp = p - tmp;
-        row[j] += tmp;
-        row[j] %= p;
+      // only the non-pivot entries have to be considered
+      // (and also the first entry)
+      row[piv] = 0;
+
+      int smallestNonPivIndex = 0;
+      while (nonPivots[smallestNonPivIndex] < piv)
+      {
+        smallestNonPivIndex++;
+      }
+
+      for (int j = smallestNonPivIndex; j < n-rows; j++)
+      {
+        unsigned ind = nonPivots[j];
+        if (matrix[i][ind] != 0)
+        {
+          unsigned long tmp = multMod (matrix[i][ind], x, p);
+          tmp = p - tmp;
+          row[ind] += tmp;
+          if (row[ind] >= p)
+          {
+            row[ind] -= p;
+          }
+        }
       }
     }
@@ -240 +274 @@
   if(piv != -1)
   {
-    // normalize and insert row into the matrix
+    // Normalize and insert row into the matrix.
+    // Then reduce upwards.
     normalizeRow (row, piv);
     for(int i = 0; i < n; i++)
@@ -247 +282 @@
     }
 
+
+    for (int i = 0; i < rows; i++)
+    {
+      unsigned x = matrix[i][piv];
+      // if the corresponding entry in the matrix is zero,
+      // there is nothing to do
+      if (x != 0)
+      {
+        for (int j = piv; j < n; j++)
+        {
+          if (row[j] != 0)
+          {
+            unsigned long tmp = multMod(row[j], x, p);
+            tmp = p - tmp;
+            matrix[i][j] += tmp;
+            if (matrix[i][j] >= p)
+            {
+              matrix[i][j] -= p;
+            }
+          }
+        }
+      }
+    }
+
     pivots[rows] = piv;
+
+    // update nonPivots
+    for (int i = 0; i < n-rows; i++)
+    {
+      if (nonPivots[i] == piv)
+      {
+        // shift everything one position to the left
+        for (int j = i; j < n-rows-1; j++)
+        {
+          nonPivots[j] = nonPivots[j+1];
+        }
+
+        break;
+      }
+    }
+
     rows++;
   }
@@ -255 +330 @@
 void NewVectorMatrix::insertMatrix (LinearDependencyMatrix & mat)
 {
-  // The matrix in LinearDependencyMatrix is already in reduced form.
-  // Thus, if the current matrix is empty, we can simply copy the other matrix.
-  if(rows == 0)
-  {
-    for(int i = 0; i < mat.rows; i++)
-    {
-      for(int j = 0; j < n; j++)
-      {
-        matrix[i][j] = mat.matrix[i][j];
-
-        rows = mat.rows;
-        for(int i = 0; i < rows; i++)
-        {
-          pivots[i] = mat.pivots[i];
-        }
-      }
-    }
-  }
-  else
-  {
-    for(int i = 0; i < mat.rows; i++)
-    {
-      insertRow (mat.matrix[i]);
-    }
+  for(int i = 0; i < mat.rows; i++)
+  {
+    insertRow (mat.matrix[i]);
   }
 }
 
 int NewVectorMatrix::findSmallestNonpivot ()
+{
+  // This method isn't very efficient, but it is called at most a few times,
+  // so efficiency is not important.
+  if(rows == n)
+    return -1;
+
+  for(int i = 0; i < n; i++)
+  {
+    bool isPivot = false;
+    for(int j = 0; j < rows; j++)
+    {
+      if(pivots[j] == i)
+      {
+        isPivot = true;
+        break;
+      }
+    }
+
+    if(!isPivot)
+    {
+      return i;
+    }
+  }
+}
+
+int NewVectorMatrix::findLargestNonpivot ()
 {
   // This method isn't very efficient, but it is called at most a few times, so efficiency is not important.
@@ -288 +368 @@
     return -1;
 
-  for(int i = 0; i < n; i++)
+  for(int i = n-1; i >= 0; i--)
   {
     bool isPivot = false;
@@ -309 +389 @@
 
 void vectorMatrixMult (unsigned long *vec, unsigned long **mat,
+                       unsigned **nonzeroIndices, unsigned *nonzeroCounts,
                        unsigned long *result, unsigned n, unsigned long p)
 {
   unsigned long tmp;
+
   for(int i = 0; i < n; i++)
   {
     result[i] = 0;
-    for(int j = 0; j < n; j++)
-    {
-      tmp = multMod (vec[j], mat[j][i], p);
+    for(int j = 0; j < nonzeroCounts[i]; j++)
+    {
+      tmp = multMod (vec[nonzeroIndices[i][j]], mat[nonzeroIndices[i][j]][i], p);
       result[i] += tmp;
-      result[i] %= p;
+      if (result[i] >= p)
+      {
+        result[i] -= p;
+      }
     }
   }
@@ -358 +443 @@
 
 
-  int i = 0;
+  // Store the indices where the matrix has non-zero entries.
+  // This has a huge impact on spares matrices.
+  unsigned* nonzeroCounts = new unsigned[n];
+  unsigned** nonzeroIndices = new unsigned*[n];
+  for (int i = 0; i < n; i++)
+  {
+    nonzeroIndices[i] = new unsigned[n];
+    nonzeroCounts[i] = 0;
+    for (int j = 0; j < n; j++)
+    {
+      if (matrix[j][i] != 0)
+      {
+        nonzeroIndices[i][nonzeroCounts[i]] = j;
+        nonzeroCounts[i]++;
+      }
+    }
+  }
+
+  int i = n-1;
 
   unsigned long *vec = new unsigned long[n];
   unsigned long *vecnew = new unsigned long[n];
 
+  unsigned loopsEven = true;
   while(i != -1)
   {
@@ -382 +486 @@
       }
 
-      vectorMatrixMult (vec, matrix, vecnew, n, p);
+      vectorMatrixMult (vec, matrix, nonzeroIndices, nonzeroCounts, vecnew, n, p);
       unsigned long *swap = vec;
       vec = vecnew;
@@ -426 +530 @@
       {
         newvectormat.insertMatrix (lindepmat);
-        i = newvectormat.findSmallestNonpivot ();
-      }
-    }
-  }
-
-  // TODO: take lcms of the different monomials!
+
+        // choose new unit vector from the front or the end, alternating
+        // for each round. If the matrix is the companion matrix for x^n,
+        // then taking vectors from the end is best. If it is the transpose,
+        // taking vectors from the front is best.
+        // This tries to take the middle way
+        if (loopsEven)
+        {
+          i = newvectormat.findSmallestNonpivot ();
+        }
+        else
+        {
+          i = newvectormat.findLargestNonpivot ();
+        }
+      }
+    }
+
+    loopsEven = !loopsEven;
+  }
+
+  for (int i = 0; i < n; i++)
+  {
+    delete[] nonzeroIndices[i];
+  }
+  delete[] nonzeroIndices;
+  delete[] nonzeroCounts;
+
 
   delete[]vecnew;
@@ -453 +578 @@
       long tmp = p - multMod (factor, q[i], p);
       a[d + i] += tmp;
-      a[d + i] %= p;
-    }
-
-    while(a[dega] == 0 && dega >= 0)
+      if (a[d + i] >= p)
+      {
+        a[d + i] -= p;
+      }
+    }
+
+    while(dega >= 0 && a[dega] == 0)
     {
       dega--;
@@ -470 +598 @@
   unsigned long *result = new unsigned long[degres + 1];
 
+  // initialize to zero
+  for (int i = 0; i <= degres; i++)
+  {
+    result[i] = 0;
+  }
+
   while(degq <= dega)
   {
     unsigned d = dega - degq;
-    long inv = modularInverse (q[degq], p);
+    unsigned long inv = modularInverse (q[degq], p);
     result[d] = multMod (a[dega], inv, p);
     for(int i = degq; i >= 0; i--)
     {
-      long tmp = p - multMod (result[d], q[i], p);
+      unsigned long tmp = p - multMod (result[d], q[i], p);
       a[d + i] += tmp;
-      a[d + i] %= p;
-    }
-
-    while(a[dega] == 0 && dega >= 0)
+      if (a[d + i] >= p)
+      {
+        a[d + i] -= p;
+      }
+    }
+
+    while(dega >= 0 && a[dega] == 0)
     {
       dega--;
@@ -515 +652 @@
     {
       result[i + j] += multMod (a[i], b[j], p);
-      result[i + j] %= p;
+      if (result[i + j] >= p)
+      {
+        result[i + j] -= p;
+      }
     }
   }

Singular/minpoly.h

@@ -1 +1 @@
 /***********************************************************************************
- * Copyright (C) 2011 Sebastian Jambor                                             *
+ * Author: Sebastian Jambor, 2011                                                  *
+ * (C) GPL (e-mail from June 6, 2012, 17:00:31 MESZ)                               *
  * sebastian@momo.math.rwth-aachen.de                                              *
  *                                                                                 *
@@ -110 +111 @@
     unsigned long **matrix;
     unsigned *pivots;
+    unsigned *nonPivots;
     unsigned rows;
 
@@ -135 +137 @@
     // If no such number exists, return -1.
     int findSmallestNonpivot();
+
+    int findLargestNonpivot();
 };
 
@@ -161 +165 @@
 unsigned long modularInverse(long long x, long long p);
 
-void vectorMatrixMult(unsigned long* vec, unsigned long **mat, unsigned long* result, unsigned n, unsigned long p);
+void vectorMatrixMult(unsigned long* vec, unsigned long **mat, unsigned **nonzeroIndices, unsigned *nonzeroCounts, unsigned long* result, unsigned n, unsigned long p);
 
 // a is a vector of length at least dega + 1, and q is a vector of length at least degq + 1,