Changeset 5b2745 in git

Timestamp:

May 27, 2011, 6:39:50 PM (13 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')

Children:

7d1c995f987c51bad98f33e83e328a16ded060ec

Parents:

db83bf13b62d944de021f52de5ae04b7ffd0390e

Message:

fix unsigned/signed int usage

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

File:

• Singular/minpoly.cc (modified) (4 diffs)

Singular/minpoly.cc

@@ -6 +6 @@
 
 #include<cmath>
+#include<Singular/mod2.h>
 //#include<iomanip>
 
@@ -12 +13 @@
 // TODO: avoid code copying, use subclassing instead!
 
-LinearDependencyMatrix::LinearDependencyMatrix(unsigned n, unsigned long p) {
-    this->n = n;
-    this->p = p;
-
-    matrix = new unsigned long*[n];
-    for (int i = 0; i < n; i++) {
-        matrix[i] = new unsigned long[2*n + 1];
-    }
-    pivots = new unsigned[n];
-    tmprow = new unsigned long[2*n + 1];
-    rows = 0;
-}
-
-LinearDependencyMatrix::~LinearDependencyMatrix() {
-    delete[] tmprow;
-    delete[] pivots;
-
-    for (int i = 0; i < n; i++) {
-        delete[] (matrix[i]);
-    }
-    delete[] matrix;
-}
-
-void LinearDependencyMatrix::resetMatrix() {
-    rows = 0;
-}
-
-int LinearDependencyMatrix::firstNonzeroEntry(unsigned long *row) {
-    for (int i = 0; i < n; i++)
-        if (row[i] != 0)
-            return i;
-
-    return -1;
-}
-
-void LinearDependencyMatrix::reduceTmpRow() {
-    for (int i = 0; i < rows; i++) {
-        unsigned piv = pivots[i];
-        unsigned long x = tmprow[piv] % p;
-        // if the corresponding entry in the row is zero, there is nothing to do
-        if (x != 0) {
-            // subtract tmprow[i] times the i-th row
-            for (int j = piv; j < n + rows + 1; j++) {
-                unsigned long tmp = matrix[i][j]*x % p;
-                tmp = p - tmp;
-                tmprow[j] += tmp;
-                // We don't normalize here, so remember to do it after all reductions.
-                // tmprow[j] %= p;
-            }
-        }
-    }
-
-    // normalize the entries of tmprow.
-    for (int i = 0; i < n + rows + 1; i++) {
-        tmprow[i] %= p;
-    }
-}
-
-
-void LinearDependencyMatrix::normalizeTmp(unsigned i) {
-    unsigned long inv = modularInverse(tmprow[i], p);
-    tmprow[i] = 1;
-    for (int j = i+1; j < 2*n + 1; j++)
-        tmprow[j] = (tmprow[j] * inv) % p;
-}
-
-bool LinearDependencyMatrix::findLinearDependency(unsigned long* newRow, unsigned long* dep) {
-    // Copy newRow to tmprow (we need to add RHS)
-    for (int i = 0; i < n; i++) {
-        tmprow[i] = newRow[i];
-        tmprow[n + i] = 0;
-    }
-    tmprow[2*n] = 0;
-    tmprow[n + rows] = 1;
-
-    reduceTmpRow();
-
-    // Is tmprow reduced to zero? Then we have found a linear dependence.
-    // Otherwise, we just add tmprow to the matrix.
-    unsigned newpivot = firstNonzeroEntry(tmprow);
-    if (newpivot == -1) {
-        for (int i = 0; i <=n; i++) {
-            dep[i] = tmprow[n+i];
-        }
-
-        return true;
-    } else {
-        normalizeTmp(newpivot);
-
-        for (int i = 0; i < 2*n + 1; i++) {
-            matrix[rows][i] = tmprow[i];
-        }
-
-        pivots[rows] = newpivot;
-        rows++;
-
-        return false;
-    }
+LinearDependencyMatrix::LinearDependencyMatrix (unsigned n, unsigned long p)
+{
+  this->n = n;
+  this->p = p;
+
+  matrix = new unsigned long *[n];
+  for(unsigned i = 0; i < n; i++)
+  {
+    matrix[i] = new unsigned long[2 * n + 1];
+  }
+  pivots = new unsigned[n];
+  tmprow = new unsigned long[2 * n + 1];
+  rows = 0;
+}
+
+LinearDependencyMatrix::~LinearDependencyMatrix ()
+{
+  delete[]tmprow;
+  delete[]pivots;
+
+  for(unsigned i = 0; i < n; i++)
+  {
+    delete[](matrix[i]);
+  }
+  delete[]matrix;
+}
+
+void LinearDependencyMatrix::resetMatrix ()
+{
+  rows = 0;
+}
+
+int LinearDependencyMatrix::firstNonzeroEntry (unsigned long *row)
+{
+  for(unsigned i = 0; i < n; i++)
+    if(row[i] != 0)
+      return i;
+
+  return -1;
+}
+
+void LinearDependencyMatrix::reduceTmpRow ()
+{
+  for(unsigned i = 0; i < rows; i++)
+  {
+    unsigned piv = pivots[i];
+    unsigned long x = tmprow[piv] % p;
+    // if the corresponding entry in the row is zero, there is nothing to do
+    if(x != 0)
+    {
+      // subtract tmprow[i] times the i-th row
+      for(unsigned j = piv; j < n + rows + 1; j++)
+      {
+        unsigned long tmp = matrix[i][j] * x % p;
+        tmp = p - tmp;
+        tmprow[j] += tmp;
+        // We don't normalize here, so remember to do it after all reductions.
+        // tmprow[j] %= p;
+      }
+    }
+  }
+
+  // normalize the entries of tmprow.
+  for(unsigned i = 0; i < n + rows + 1; i++)
+  {
+    tmprow[i] %= p;
+  }
+}
+
+
+void LinearDependencyMatrix::normalizeTmp (unsigned i)
+{
+  unsigned long inv = modularInverse (tmprow[i], p);
+  tmprow[i] = 1;
+  for(unsigned j = i + 1; j < 2 * n + 1; j++)
+    tmprow[j] = (tmprow[j] * inv) % p;
+}
+
+bool LinearDependencyMatrix::findLinearDependency (unsigned long *newRow,
+                                                   unsigned long *dep)
+{
+  // Copy newRow to tmprow (we need to add RHS)
+  for(unsigned i = 0; i < n; i++)
+  {
+    tmprow[i] = newRow[i];
+    tmprow[n + i] = 0;
+  }
+  tmprow[2 * n] = 0;
+  tmprow[n + rows] = 1;
+
+  reduceTmpRow ();
+
+  // Is tmprow reduced to zero? Then we have found a linear dependence.
+  // Otherwise, we just add tmprow to the matrix.
+  int newpivot = firstNonzeroEntry (tmprow);
+  if(newpivot == -1)
+  {
+    for(unsigned i = 0; i <= n; i++)
+    {
+      dep[i] = tmprow[n + i];
+    }
+
+    return true;
+  }
+  else
+  {
+    normalizeTmp (newpivot);
+
+    for(unsigned i = 0; i < 2 * n + 1; i++)
+    {
+      matrix[rows][i] = tmprow[i];
+    }
+
+    pivots[rows] = newpivot;
+    rows++;
+
+    return false;
+  }
 }
 
 #if 0
-std::ostream& operator<<(std::ostream& out, const LinearDependencyMatrix& mat) {
-    int width = ((int)log10(mat.p)) + 1;
-
-    for (int i = 0; i < mat.rows; i++) {
-        for (int j = 0; j < mat.n; j++) {
-            out << std::setw(width) << mat.matrix[i][j] << " ";
-        }
-        out << "| ";
-        for (int j = mat.n; j <= 2*mat.n; j++) {
-            out << std::setw(width) << mat.matrix[i][j] << " ";
-        }
-        out << std::endl;
-    }
-    out << "-----" << std::endl;
-    for (int j = 0; j < mat.n; j++) {
-        out << std::setw(width) << mat.tmprow[j] << " ";
+std::ostream & operator<< (std::ostream & out,
+                           const LinearDependencyMatrix & mat)
+{
+  int width = ((int) log10 (mat.p)) + 1;
+
+  for(int i = 0; i < mat.rows; i++)
+  {
+    for(int j = 0; j < mat.n; j++)
+    {
+      out << std::setw (width) << mat.matrix[i][j] << " ";
     }
     out << "| ";
-    for (int j = mat.n; j <= 2*mat.n; j++) {
-        out << std::setw(width) << mat.tmprow[j] << " ";
+    for(int j = mat.n; j <= 2 * mat.n; j++)
+    {
+      out << std::setw (width) << mat.matrix[i][j] << " ";
     }
     out << std::endl;
-
-    return out;
+  }
+  out << "-----" << std::endl;
+  for(int j = 0; j < mat.n; j++)
+  {
+    out << std::setw (width) << mat.tmprow[j] << " ";
+  }
+  out << "| ";
+  for(int j = mat.n; j <= 2 * mat.n; j++)
+  {
+    out << std::setw (width) << mat.tmprow[j] << " ";
+  }
+  out << std::endl;
+
+  return out;
 }
 #endif
 
 
-NewVectorMatrix::NewVectorMatrix(unsigned n, unsigned long p) {
-    this->n = n;
-    this->p = p;
-
-    matrix = new unsigned long*[n];
-    for (int i = 0; i < n; i++) {
-        matrix[i] = new unsigned long[n];
-    }
-
-    pivots = new unsigned[n];
-    rows = 0;
-}
-
-NewVectorMatrix::~NewVectorMatrix() {
-    delete pivots;
-
-    for (int i = 0; i < n; i++) {
-        delete[] matrix[i];
-    }
-    delete matrix;
-}
-
-int NewVectorMatrix::firstNonzeroEntry(unsigned long *row) {
-    for (int i = 0; i < n; i++)
-        if (row[i] != 0)
-            return i;
-
-    return -1;
+NewVectorMatrix::NewVectorMatrix (unsigned n, unsigned long p)
+{
+  this->n = n;
+  this->p = p;
+
+  matrix = new unsigned long *[n];
+  for(unsigned i = 0; i < n; i++)
+  {
+    matrix[i] = new unsigned long[n];
+  }
+
+  pivots = new unsigned[n];
+  rows = 0;
+}
+
+NewVectorMatrix::~NewVectorMatrix ()
+{
+  delete pivots;
+
+  for(unsigned i = 0; i < n; i++)
+  {
+    delete[]matrix[i];
+  }
+  delete matrix;
+}
+
+int NewVectorMatrix::firstNonzeroEntry (unsigned long *row)
+{
+  for(unsigned i = 0; i < n; i++)
+    if(row[i] != 0)
+      return i;
+
+  return -1;
 }
 
@@ -182 +215 @@
 //}
 //
-void NewVectorMatrix::normalizeRow(unsigned long *row, unsigned i) {
-    unsigned long inv = modularInverse(row[i], p);
-    row[i] = 1;
-
-    for (int j = i+1; j < n; j++) {
-        row[j] = (row[j] * inv) % p;
-    }
-}
-
-void NewVectorMatrix::insertRow(unsigned long* row) {
-    for (int i = 0; i < rows; i++) {
-        unsigned piv = pivots[i];
-        unsigned long x = row[piv] % p;
-        // if the corresponding entry in the row is zero, there is nothing to do
-        if (x != 0) {
-            // subtract row[i] times the i-th row
-            for (int j = piv; j < n; j++) {
-                unsigned long tmp = matrix[i][j]*x % p;
-                tmp = p - tmp;
-                row[j] += tmp;
-                // We don't normalize here, so remember to do it after all reductions.
-                // row[j] %= p;
-            }
+void NewVectorMatrix::normalizeRow (unsigned long *row, unsigned i)
+{
+  unsigned long inv = modularInverse (row[i], p);
+  row[i] = 1;
+
+  for(unsigned j = i + 1; j < n; j++)
+  {
+    row[j] = (row[j] * inv) % p;
+  }
+}
+
+void NewVectorMatrix::insertRow (unsigned long *row)
+{
+  unsigned i;
+  for(i = 0; i < rows; i++)
+  {
+    unsigned piv = pivots[i];
+    unsigned long x = row[piv] % p;
+    // if the corresponding entry in the row is zero, there is nothing to do
+    if(x != 0)
+    {
+      // subtract row[i] times the i-th row
+      for(int j = piv; j < (int)n; j++)
+      {
+        unsigned long tmp = matrix[i][j] * x % p;
+        tmp = p - tmp;
+        row[j] += tmp;
+        // We don't normalize here, so remember to do it after all reductions.
+        // row[j] %= p;
+      }
+    }
+  }
+
+  // normalize the entries of row.
+  for(i = 0; i < n + rows + 1; i++)
+  {
+    row[i] %= p;
+  }
+
+  int piv = firstNonzeroEntry (row);
+
+  if(piv != -1)
+  {
+    // normalize and insert row into the matrix
+    normalizeRow (row, piv);
+    for(i = 0; i < n; i++)
+    {
+      matrix[rows][i] = row[i];
+    }
+
+    pivots[rows] = piv;
+    rows++;
+  }
+}
+
+
+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(unsigned i = 0; i < mat.rows; i++)
+    {
+      for(unsigned j = 0; j < n; j++)
+      {
+        matrix[i][j] = mat.matrix[i][j];
+
+        rows = mat.rows;
+        for(unsigned ii = 0; ii < rows; ii++)
+        {
+          pivots[ii] = mat.pivots[ii];
         }
-    }
-
-    // normalize the entries of row.
-    for (int i = 0; i < n + rows + 1; i++) {
-        row[i] %= p;
-    }
-
-    unsigned piv = firstNonzeroEntry(row);
-
-    if (piv != -1) {
-        // normalize and insert row into the matrix
-        normalizeRow(row, piv);
-        for (int i = 0; i < n; i++) {
-            matrix[rows][i] = row[i];
-        }
-
-        pivots[rows] = piv;
-        rows++;
-    }
-}
-
-
-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]);
-        }
-    }
-}
-
-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;
-        }
-    }
-}
-
-
-void vectorMatrixMult(unsigned long* vec, unsigned long **mat, 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 = vec[j]*mat[j][i] % p;
-            result[i] += tmp;
-        }
-        // We can afford to reduce mod p only after all additions, since p < 2^31, but an unsigned long can store 2^64.
-        // Thus the only possibility for an overflow would occurr for matrices with about 2^31 rows.
-        result[i] %= p;
-    }
+      }
+    }
+  }
+  else
+  {
+    for(unsigned 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(unsigned i = 0; i < n; i++)
+  {
+    bool isPivot = false;
+    for(unsigned j = 0; j < rows; j++)
+    {
+      if(pivots[j] == i)
+      {
+        isPivot = true;
+        break;
+      }
+    }
+
+    if(!isPivot)
+    {
+      return i;
+    }
+  }
+  assume(0);
+  return -1; // to make the compiler happy
+}
+
+
+void vectorMatrixMult (unsigned long *vec, unsigned long **mat,
+                       unsigned long *result, unsigned n, unsigned long p)
+{
+  unsigned long tmp;
+  for(unsigned i = 0; i < n; i++)
+  {
+    result[i] = 0;
+    for(unsigned j = 0; j < n; j++)
+    {
+      tmp = vec[j] * mat[j][i] % p;
+      result[i] += tmp;
+    }
+    // We can afford to reduce mod p only after all additions, since p < 2^31, but an unsigned long can store 2^64.
+    // Thus the only possibility for an overflow would occurr for matrices with about 2^31 rows.
+    result[i] %= p;
+  }
 }
 
 #if 0
-void printVec(unsigned long* vec, int n) {
-    for (int i = 0; i < n; i++) {
-        std::cout << vec[i] << ", ";
-    }
-
-    std::cout << std::endl;
+void printVec (unsigned long *vec, int n)
+{
+  for(int i = 0; i < n; i++)
+  {
+    std::cout << vec[i] << ", ";
+  }
+
+  std::cout << std::endl;
 }
 #endif
 
 
-unsigned long* computeMinimalPolynomial(unsigned long** matrix, unsigned n, unsigned long p) {
-    LinearDependencyMatrix lindepmat(n, p);
-    NewVectorMatrix newvectormat(n, p);
-
-    unsigned long* result = new unsigned long[n+1];
-    unsigned long* mpvec = new unsigned long[n+1];
-    unsigned long* tmp = new unsigned long[n+1];
-
-    // initialize result = 1
-    for (int i = 0; i <= n; i++) {
-        result[i] = 0;
-    }
-    result[0] = 1;
-
-    int degresult = 0;
-
-
-    int i = 0;
-
-    unsigned long* vec = new unsigned long[n];
-    unsigned long* vecnew = new unsigned long[n];
-
-    while (i != -1) {
-        for (int j = 0; j < n; j++) {
-            vec[j] = 0;
-        }
-        vec[i] = 1;
-
-        lindepmat.resetMatrix();
-
-        while (true) {
-            bool ld = lindepmat.findLinearDependency(vec, mpvec);
-
-            if (ld) {
-                break;
-            }
-
-            vectorMatrixMult(vec, matrix, vecnew, n, p);
-            unsigned long* swap = vec;
-            vec = vecnew;
-            vecnew = swap;
-        }
-
-
-        unsigned degmpvec = n;
-        while (mpvec[degmpvec] == 0) {
-            degmpvec--;
-        }
-
-        // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree,
-        // then we are done.
-        if (degmpvec == n) {
-            unsigned long* swap = result;
-            result = mpvec;
-            mpvec = swap;
-            i = -1;
-        } else {
-            // otherwise, we have to compute the lcm of mpvec and prior result
-
-            // tmp = zeropol
-            for (int j = 0; j <= n; j++) {
-                tmp[j] = 0;
-            }
-            degresult  = lcm(tmp, result, mpvec, p, degresult, degmpvec);
-            unsigned long* swap = result;
-            result = tmp;
-            tmp = swap;
-
-            if (degresult == n) {
-                // minimal polynomial has highest possible degree, so we are done.
-                i = -1;
-            } else {
-                newvectormat.insertMatrix(lindepmat);
-                i = newvectormat.findSmallestNonpivot();
-            }
-        }
-    }
-
-    // TODO: take lcms of the different monomials!
-
-    delete[] vecnew;
-    delete[] vec;
-    delete[] tmp;
-    delete[] mpvec;
-
-    return result;
-}
-
-
-void rem(unsigned long* a, unsigned long* q, unsigned long p, int & dega, int degq) {
-    while (degq <= dega) {
-        unsigned d = dega - degq;
-        long factor = a[dega]*modularInverse(q[degq], p) % p;
-        for (int i = degq; i >= 0; i--) {
-            long tmp = p - (factor*q[i] % p);
-            a[d + i] += tmp;
-            a[d + i] %= p;
-        }
-
-        while (a[dega] == 0 && dega >= 0) {
-            dega--;
-        }
-    }
-}
-
-
-void quo(unsigned long* a, unsigned long* q, unsigned long p, int & dega, int degq) {
-    unsigned degres = dega - degq;
-    unsigned long* result = new unsigned long[degres + 1];
-
-    while (degq <= dega) {
-        unsigned d = dega - degq;
-        long inv = modularInverse(q[degq], p);
-        result[d] = a[dega]*inv % p;
-        for (int i = degq; i >= 0; i--) {
-            long tmp = p - (result[d]*q[i] % p);
-            a[d + i] += tmp;
-            a[d + i] %= p;
-        }
-
-        while (a[dega] == 0 && dega >= 0) {
-            dega--;
-        }
-    }
-
-    // copy result into a
-    for (int i = 0; i <= degres; i++) {
-        a[i] = result[i];
-    }
-    // set all following entries of a to zero
-    for (int i = degres + 1; i <= degq + degres; i++) {
-        a[i] = 0;
-    }
-
-    dega = degres;
-
-    delete[] result;
-}
-
-
-void mult(unsigned long* result, unsigned long* a, unsigned long* b, unsigned long p, int dega, int degb) {
-    // NOTE: we assume that every entry in result is preinitialized to zero!
-
-    for (int i = 0; i <= dega; i++) {
-        for (int j = 0; j <= degb; j++) {
-            result[i + j] += (a[i]*b[j] % p);
-            result[i + j] %= p;
-        }
-    }
-}
-
-
-int gcd(unsigned long* g, unsigned long* a, unsigned long* b, unsigned long p, int dega, int degb) {
-    unsigned long* tmp1 = new unsigned long[dega + 1];
-    unsigned long* tmp2 = new unsigned long[degb + 1];
-    for (int i = 0; i <= dega; i++) {
-        tmp1[i] = a[i];
-    }
-    for (int i = 0; i <= degb; i++) {
-        tmp2[i] = b[i];
-    }
-    int degtmp1 = dega;
-    int degtmp2 = degb;
-
-    unsigned long* swappol;
-    int swapdeg;
-
-    while (degtmp2 >= 0) {
-        rem(tmp1, tmp2, p, degtmp1, degtmp2);
-        swappol = tmp1;
-        tmp1 = tmp2;
-        tmp2 = swappol;
-
-        swapdeg = degtmp1;
-        degtmp1 = degtmp2;
-        degtmp2 = swapdeg;
-    }
-
-    for (int i = 0; i <= degtmp1; i++) {
-        g[i] = tmp1[i];
-    }
-
-    delete[] tmp1;
-    delete[] tmp2;
-
-    return degtmp1;
-}
-
-
-int lcm(unsigned long* l, unsigned long* a, unsigned long* b, unsigned long p, int dega, int degb) {
-    unsigned long* g = new unsigned long[dega + 1];
-    // initialize entries of g to zero
-    for (int i = 0; i <= dega; i++) {
-        g[i] = 0;
-    }
-
-    int degg = gcd(g, a, b, p, dega, degb);
-
-    if (degg > 0) {
-        // non-trivial gcd, so compute a = (a/g)
-        quo(a, g, p, dega, degg);
-    }
-    mult(l, a, b, p, dega, degb);
-
-    // normalize
-    if (l[dega + degb + 1] != 1) {
-        unsigned long inv = modularInverse(l[dega + degb], p);
-        for (int i = 0; i <= dega + degb; i++) {
-            l[i] *= inv;
-            l[i] %= p;
-        }
-    }
-
-    return dega + degb;
+unsigned long *computeMinimalPolynomial (unsigned long **matrix, unsigned n,
+                                         unsigned long p)
+{
+  LinearDependencyMatrix lindepmat (n, p);
+  NewVectorMatrix newvectormat (n, p);
+
+  unsigned long *result = new unsigned long[n + 1];
+  unsigned long *mpvec = new unsigned long[n + 1];
+  unsigned long *tmp = new unsigned long[n + 1];
+
+  // initialize result = 1
+  for(unsigned i = 0; i <= n; i++)
+  {
+    result[i] = 0;
+  }
+  result[0] = 1;
+
+  int degresult = 0;
+
+
+  int i = 0;
+
+  unsigned long *vec = new unsigned long[n];
+  unsigned long *vecnew = new unsigned long[n];
+
+  while(i != -1)
+  {
+    for(unsigned j = 0; j < n; j++)
+    {
+      vec[j] = 0;
+    }
+    vec[i] = 1;
+
+    lindepmat.resetMatrix ();
+
+    while(true)
+    {
+      bool ld = lindepmat.findLinearDependency (vec, mpvec);
+
+      if(ld)
+      {
+        break;
+      }
+
+      vectorMatrixMult (vec, matrix, vecnew, n, p);
+      unsigned long *swap = vec;
+      vec = vecnew;
+      vecnew = swap;
+    }
+
+
+    unsigned degmpvec = n;
+    while(mpvec[degmpvec] == 0)
+    {
+      degmpvec--;
+    }
+
+    // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree,
+    // then we are done.
+    if(degmpvec == n)
+    {
+      unsigned long *swap = result;
+      result = mpvec;
+      mpvec = swap;
+      i = -1;
+    }
+    else
+    {
+      // otherwise, we have to compute the lcm of mpvec and prior result
+
+      // tmp = zeropol
+      for(unsigned j = 0; j <= n; j++)
+      {
+        tmp[j] = 0;
+      }
+      degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec);
+      unsigned long *swap = result;
+      result = tmp;
+      tmp = swap;
+
+      if(degresult == n)
+      {
+        // minimal polynomial has highest possible degree, so we are done.
+        i = -1;
+      }
+      else
+      {
+        newvectormat.insertMatrix (lindepmat);
+        i = newvectormat.findSmallestNonpivot ();
+      }
+    }
+  }
+
+  // TODO: take lcms of the different monomials!
+
+  delete[]vecnew;
+  delete[]vec;
+  delete[]tmp;
+  delete[]mpvec;
+
+  return result;
+}
+
+
+void rem (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
+          int degq)
+{
+  while(degq <= dega)
+  {
+    unsigned d = dega - degq;
+    long factor = a[dega] * modularInverse (q[degq], p) % p;
+    for(int i = degq; i >= 0; i--)
+    {
+      long tmp = p - (factor * q[i] % p);
+      a[d + i] += tmp;
+      a[d + i] %= p;
+    }
+
+    while(a[dega] == 0 && dega >= 0)
+    {
+      dega--;
+    }
+  }
+}
+
+
+void quo (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
+          int degq)
+{
+  unsigned degres = dega - degq;
+  unsigned long *result = new unsigned long[degres + 1];
+  unsigned i;
+
+  while(degq <= dega)
+  {
+    unsigned d = dega - degq;
+    long inv = modularInverse (q[degq], p);
+    result[d] = a[dega] * inv % p;
+    for(i = degq; i >= 0; i--)
+    {
+      long tmp = p - (result[d] * q[i] % p);
+      a[d + i] += tmp;
+      a[d + i] %= p;
+    }
+
+    while(a[dega] == 0 && dega >= 0)
+    {
+      dega--;
+    }
+  }
+
+  // copy result into a
+  for(i = 0; i <= degres; i++)
+  {
+    a[i] = result[i];
+  }
+  // set all following entries of a to zero
+  for(i = degres + 1; i <= degq + degres; i++)
+  {
+    a[i] = 0;
+  }
+
+  dega = degres;
+
+  delete[]result;
+}
+
+
+void mult (unsigned long *result, unsigned long *a, unsigned long *b,
+           unsigned long p, int dega, int degb)
+{
+  // NOTE: we assume that every entry in result is preinitialized to zero!
+
+  for(int i = 0; i <= dega; i++)
+  {
+    for(int j = 0; j <= degb; j++)
+    {
+      result[i + j] += (a[i] * b[j] % p);
+      result[i + j] %= p;
+    }
+  }
+}
+
+
+int gcd (unsigned long *g, unsigned long *a, unsigned long *b,
+         unsigned long p, int dega, int degb)
+{
+  unsigned long *tmp1 = new unsigned long[dega + 1];
+  unsigned long *tmp2 = new unsigned long[degb + 1];
+  for(int i = 0; i <= dega; i++)
+  {
+    tmp1[i] = a[i];
+  }
+  for(int i = 0; i <= degb; i++)
+  {
+    tmp2[i] = b[i];
+  }
+  int degtmp1 = dega;
+  int degtmp2 = degb;
+
+  unsigned long *swappol;
+  int swapdeg;
+
+  while(degtmp2 >= 0)
+  {
+    rem (tmp1, tmp2, p, degtmp1, degtmp2);
+    swappol = tmp1;
+    tmp1 = tmp2;
+    tmp2 = swappol;
+
+    swapdeg = degtmp1;
+    degtmp1 = degtmp2;
+    degtmp2 = swapdeg;
+  }
+
+  for(int i = 0; i <= degtmp1; i++)
+  {
+    g[i] = tmp1[i];
+  }
+
+  delete[]tmp1;
+  delete[]tmp2;
+
+  return degtmp1;
+}
+
+
+int lcm (unsigned long *l, unsigned long *a, unsigned long *b,
+         unsigned long p, int dega, int degb)
+{
+  unsigned long *g = new unsigned long[dega + 1];
+  // initialize entries of g to zero
+  for(int i = 0; i <= dega; i++)
+  {
+    g[i] = 0;
+  }
+
+  int degg = gcd (g, a, b, p, dega, degb);
+
+  if(degg > 0)
+  {
+    // non-trivial gcd, so compute a = (a/g)
+    quo (a, g, p, dega, degg);
+  }
+  mult (l, a, b, p, dega, degb);
+
+  // normalize
+  if(l[dega + degb + 1] != 1)
+  {
+    unsigned long inv = modularInverse (l[dega + degb], p);
+    for(int i = 0; i <= dega + degb; i++)
+    {
+      l[i] *= inv;
+      l[i] %= p;
+    }
+  }
+
+  return dega + degb;
 }
 
 
 // compute x^(-1) mod p
-long modularInverse(long x, long p) {
-    long u1 = 1;
-    long u2 = 0;
-    long u3 = x;
-    long v1 = 0;
-    long v2 = 1;
-    long v3 = p;
-
-    long q, t1, t2, t3;
-    while (v3 != 0) {
-        q = u3 / v3;
-        t1 = u1 - q * v1;
-        t2 = u2 - q * v2;
-        t3 = u3 - q * v3;
-        u1 = v1;
-        u2 = v2;
-        u3 = v3;
-        v1 = t1;
-        v2 = t2;
-        v3 = t3;
-    }
-
-    if (u1 < 0) {
-        u1 += p;
-    }
-
-    return u1;
-}
+long modularInverse (long x, long p)
+{
+  long u1 = 1;
+  long u2 = 0;
+  long u3 = x;
+  long v1 = 0;
+  long v2 = 1;
+  long v3 = p;
+
+  long q, t1, t2, t3;
+  while(v3 != 0)
+  {
+    q = u3 / v3;
+    t1 = u1 - q * v1;
+    t2 = u2 - q * v2;
+    t3 = u3 - q * v3;
+    u1 = v1;
+    u2 = v2;
+    u3 = v3;
+    v1 = t1;
+    v2 = t2;
+    v3 = t3;
+  }
+
+  if(u1 < 0)
+  {
+    u1 += p;
+  }
+
+  return u1;
+}
+
 /*
 int main() {
@@ -744 +856 @@
 }
 */
-