source: git/Singular/MinorProcessor.cc @ b08f6fe

#include "mod2.h"
#include "structs.h"
#include "polys.h"
#include <MinorProcessor.h>
#include "febase.h"
#include "kstd1.h"
#include "kbuckets.h"

void MinorProcessor::print() const {
    PrintS(this->toString().c_str());
}

int MinorProcessor::getBestLine (const int k, const MinorKey& mk) const {
    // This method identifies the row or column with the most zeros.
    // The returned index (bestIndex) is absolute within the pre-defined matrix.
    // If some row has the most zeros, then the absolute (0-based) row index is returned.
    // If, contrariwise, some column has the most zeros, then -1 minus the absolute (0-based) column index is returned.
    int numberOfZeros = 0;
    int bestIndex = 100000;    // We start with an invalid row/column index.
    int maxNumberOfZeros = -1; // We update this variable whenever we find a new so-far optimal row or column.
    for (int r = 0; r < k; r++) {
        // iterate through all k rows of the momentary minor
        int absoluteR = mk.getAbsoluteRowIndex(r);
        numberOfZeros = 0;
        for (int c = 0; c < k; c++) {
            int absoluteC = mk.getAbsoluteColumnIndex(c);
            if (isEntryZero(absoluteR, absoluteC)) numberOfZeros++;
        }
        if (numberOfZeros > maxNumberOfZeros) {
            // We found a new best line which is a row.
            bestIndex = absoluteR;
            maxNumberOfZeros = numberOfZeros;
        }
    };
    for (int c = 0; c < k; c++) {
        int absoluteC = mk.getAbsoluteColumnIndex(c);
        numberOfZeros = 0;
        for (int r = 0; r < k; r++) {
            int absoluteR = mk.getAbsoluteRowIndex(r);
            if (isEntryZero(absoluteR, absoluteC)) numberOfZeros++;
        }
        if (numberOfZeros > maxNumberOfZeros) {
            // We found a new best line which is a column. So we transform the return value.
            // Note that we can easily get back absoluteC from bestLine: absoluteC = - 1 - bestLine.
            bestIndex = - absoluteC - 1;
            maxNumberOfZeros = numberOfZeros;
        }
    };
    return bestIndex;
}

void MinorProcessor::setMinorSize(const int minorSize) {
    _minorSize = minorSize;
    _minor.reset();
}

bool MinorProcessor::hasNextMinor() {
    return setNextKeys(_minorSize);
}

void MinorProcessor::getCurrentRowIndices(int* const target) const {
    return _minor.getAbsoluteRowIndices(target);
}

void MinorProcessor::getCurrentColumnIndices(int* const target) const {
    return _minor.getAbsoluteColumnIndices(target);
}

void MinorProcessor::defineSubMatrix(const int numberOfRows, const int* rowIndices,
                                     const int numberOfColumns, const int* columnIndices) {
    // The method assumes ascending row and column indices in the two argument arrays.
    // These indices are understood to be zero-based.
    // The method will set the two arrays of ints in _container.
    // Example: The indices 0, 2, 3, 7 will be converted to an array with one int
    // representing the binary number 10001101 (check bits from right to left).

    _containerRows = numberOfRows;
    int highestRowIndex = rowIndices[numberOfRows - 1];
    int rowBlockCount = (highestRowIndex / 32) + 1;
    unsigned int rowBlocks[rowBlockCount];
    for (int i = 0; i < rowBlockCount; i++) rowBlocks[i] = 0;
    for (int i = 0; i < numberOfRows; i++) {
        int blockIndex = rowIndices[i] / 32;
        int offset = rowIndices[i] % 32;
        rowBlocks[blockIndex] += (1 << offset);
    }

    _containerColumns = numberOfColumns;
    int highestColumnIndex = columnIndices[numberOfColumns - 1];
    int columnBlockCount = (highestColumnIndex / 32) + 1;
    unsigned int columnBlocks[columnBlockCount];
    for (int i = 0; i < columnBlockCount; i++) columnBlocks[i] = 0;
    for (int i = 0; i < numberOfColumns; i++) {
        int blockIndex = columnIndices[i] / 32;
        int offset = columnIndices[i] % 32;
        columnBlocks[blockIndex] += (1 << offset);
    }

    _container.set(rowBlockCount, rowBlocks, columnBlockCount, columnBlocks);
}

bool MinorProcessor::setNextKeys(const int k) {
    // This method moves _minor to the next valid kxk-minor within _container.
    // It returns true iff this is successful, i.e. iff _minor did not already encode the final kxk-minor.
    if (_minor.compare(MinorKey(0, 0, 0, 0)) == 0) {
        // This means that we haven't started yet. Thus, we are about to compute the first kxk-minor.
        _minor.selectFirstRows(k, _container);
        _minor.selectFirstColumns(k, _container);
        return true;
    }
    else if (_minor.selectNextColumns(k, _container)) {
        // Here we were able to pick a next subset of columns (within the same subset of rows).
        return true;
    }
    else if (_minor.selectNextRows(k, _container)) {
        // Here we were not able to pick a next subset of columns (within the same subset of rows).
        // But we could pick a next subset of rows. We must hence reset the subset of columns:
        _minor.selectFirstColumns(k, _container);
        return true;
    }
    else {
        // We were neither able to pick a next subset of columns nor of rows.
        // I.e., we have iterated through all sensible choices of subsets of rows and columns.
        return false;
    }
}

bool MinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
{
  assume(false);
  return false;
}

string MinorProcessor::toString () const
{
  assume(false);
  return "";
}

int MinorProcessor::IOverJ(const int i, const int j) {
    // This is a non-recursive implementation.
    assert(i >= 0 && j >= 0 && i >= j);
    if (j == 0 || i == j) return 1;
    int result = 1;
    for (int k = i - j + 1; k <= i; k++) result *= k;
    // Here, we have result = (i - j + 1) * ... * i.
    for (int k = 2; k <= j; k++) result /= k;
    // Here, we have result = (i - j + 1) * ... * i / 1 / 2 ... / j = i! / j! / (i - j)!.
    return result;
}

int MinorProcessor::Faculty(const int i) {
    // This is a non-recursive implementation.
    assert(i >= 0);
    int result = 1;
    for (int j = 1; j <= i; j++) result *= j;
    // Here, we have result = 1 * 2 * ... * i = i!
    return result;
}

int MinorProcessor::NumberOfRetrievals (const int rows, const int columns, const int containerMinorSize,
                        const int minorSize, const bool multipleMinors) {
    // This method computes the number of potential retrievals of a single minor when computing
    // all minors of a given size within a matrix of given size.
    int result = 0;
    if (multipleMinors) {
        // Here, we would like to compute all minors of size
        // containerMinorSize x containerMinorSize in a matrix of size rows x columns.
        // Then, we need to retrieve any minor of size minorSize x minorSize exactly
        // n times, where n is as follows:
        result = IOverJ(rows - minorSize, containerMinorSize - minorSize)
               * IOverJ(columns - minorSize, containerMinorSize - minorSize)
               * Faculty(containerMinorSize - minorSize);
    }
    else {
        // Here, we would like to compute just one minor of size
        // containerMinorSize x containerMinorSize. Then, we need to retrieve
        // any minor of size minorSize x minorSize exactly
        // (containerMinorSize - minorSize)! times:
        result = Faculty(containerMinorSize - minorSize);
    }
    return result;
}

MinorProcessor::MinorProcessor () {
    _container = MinorKey(0, 0, 0, 0);
    _minor = MinorKey(0, 0, 0, 0);
    _containerRows = 0;
    _containerColumns = 0;
    _minorSize = 0;
    _rows = 0;
    _columns = 0;
}

IntMinorProcessor::IntMinorProcessor () {
    _intMatrix = 0;
}

string IntMinorProcessor::toString () const {
    char h[32];
    string t = "";
    string s = "IntMinorProcessor:";
    s += "\n   matrix: ";
    sprintf(h, "%d", _rows); s += h;
    s += " x ";
    sprintf(h, "%d", _columns); s += h;
    for (int r = 0; r < _rows; r++) {
        s += "\n      ";
        for (int c = 0; c < _columns; c++) {
            sprintf(h, "%d", getEntry(r, c)); t = h;
            for (int k = 0; k < int(4 - strlen(h)); k++) s += " ";
            s += t;
        }
    }
    int myIndexArray[500];
    s += "\n   considered submatrix has row indices: ";
    _container.getAbsoluteRowIndices(myIndexArray);
    for (int k = 0; k < _containerRows; k++) {
        if (k != 0) s += ", ";
        sprintf(h, "%d", myIndexArray[k]); s += h;
    }
    s += " (first row of matrix has index 0)";
    s += "\n   considered submatrix has column indices: ";
    _container.getAbsoluteColumnIndices(myIndexArray);
    for (int k = 0; k < _containerColumns; k++) {
        if (k != 0) s += ", ";
        sprintf(h, "%d", myIndexArray[k]); s += h;
    }
    s += " (first column of matrix has index 0)";
    s += "\n   size of considered minor(s): ";
    sprintf(h, "%d", _minorSize); s += h;
    s += "x";
    s += h;
    return s;
}

IntMinorProcessor::~IntMinorProcessor() {
    // free memory of _intMatrix
    delete [] _intMatrix; _intMatrix = 0;
}

bool IntMinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
{
  return getEntry(absoluteRowIndex, absoluteColumnIndex) == 0;
}

void IntMinorProcessor::defineMatrix (const int numberOfRows, const int numberOfColumns, const int* matrix) {
    // free memory of _intMatrix
    delete [] _intMatrix; _intMatrix = 0;

    _rows = numberOfRows;
    _columns = numberOfColumns;

    // allocate memory for new entries in _intMatrix
    int n = _rows * _columns;
    _intMatrix = new int[n];

    // copying values from one-dimensional method parameter "matrix"
    for (int i = 0; i < n; i++)
      _intMatrix[i] = matrix[i];
}

IntMinorValue IntMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
                                          Cache<MinorKey, IntMinorValue>& c, const int characteristic, const ideal& iSB) {
    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
    _minorSize = dimension;
    // call a helper method which recursively computes the minor using the cache c:
    return getMinorPrivateLaplace(dimension, _container, false, c, characteristic, iSB);
}

IntMinorValue IntMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
                                          const int characteristic, const ideal& iSB, const char* algorithm) {
    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
    _minorSize = dimension;

    // call a helper method which computes the minor (without using a cache):
    if (strcmp(algorithm, "Laplace") == 0)
      return getMinorPrivateLaplace(_minorSize, _container, characteristic, iSB);
    else if (strcmp(algorithm, "Bareiss") == 0)
      return getMinorPrivateBareiss(_minorSize, _container, characteristic, iSB);
    else assume(false);
}

IntMinorValue IntMinorProcessor::getNextMinor(const int characteristic, const ideal& iSB, const char* algorithm) {
    // call a helper method which computes the minor (without using a cache):
    if (strcmp(algorithm, "Laplace") == 0)
      return getMinorPrivateLaplace(_minorSize, _minor, characteristic, iSB);
    else if (strcmp(algorithm, "Bareiss") == 0)
      return getMinorPrivateBareiss(_minorSize, _minor, characteristic, iSB);
    else assume(false);
}

IntMinorValue IntMinorProcessor::getNextMinor(Cache<MinorKey, IntMinorValue>& c, const int characteristic, const ideal& iSB) {
    // computation with cache
    return getMinorPrivateLaplace(_minorSize, _minor, true, c, characteristic, iSB);
}

/* compute the reduction of an integer i modulo an ideal which captures a std basis */
int getReduction (const int i, const ideal& iSB)
{
  if (i == 0) return 0;
  poly f = pISet(i);
  poly g = kNF(iSB, currRing->qideal, f);
  int result = 0;
  if (g != NULL) result = n_Int(pGetCoeff(g), currRing);
  pDelete(&f);
  pDelete(&g);
  return result;
}

IntMinorValue IntMinorProcessor::getMinorPrivateLaplace(const int k, const MinorKey& mk, const int characteristic, const ideal& iSB) {
    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
    if (k == 1) {
        int e = getEntry(mk.getAbsoluteRowIndex(0), mk.getAbsoluteColumnIndex(0));
        if (characteristic != 0) e = e % characteristic;
        if (iSB != 0) e = getReduction(e, iSB);
        return IntMinorValue(e, 0, 0, 0, 0, -1, -1); // "-1" is to signal that any statistics about the
                                                     // number of retrievals does not make sense, as we do not use a cache.
    }
    else {
        // Here, the minor must be 2x2 or larger.
        int b = getBestLine(k, mk);                          // row or column with most zeros
        int result = 0;                                      // This will contain the value of the minor.
        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
                                                             // ..."a*" for accumulated operation counters
        if (b >= 0) {
            // This means that the best line is the row with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorRowKey:
            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
            for (int c = 0; c < k; c++) {
                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
                if (getEntry(b, absoluteC) != 0) { // Only then do we have to consider this sub-determinante.
                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is mk with row b and column absoluteC omitted.
                    IntMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, characteristic, iSB);  // recursive call
                    m += mv.getMultiplications();
                    s += mv.getAdditions();
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    result += sign * mv.getResult() * getEntry(b, absoluteC); // adding sub-determinante times matrix entry
                                                                              // times appropriate sign
                    if (characteristic != 0) result = result % characteristic;
                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        else {
            b = - b - 1;
            // This means that the best line is the column with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorColumnKey:
            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
            for (int r = 0; r < k; r++) {
                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
                if (getEntry(absoluteR, b) != 0) { // Only then do we have to consider this sub-determinante.
                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is mk with row absoluteR and column b omitted.
                    IntMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, characteristic, iSB);  // recursive call
                    m += mv.getMultiplications();
                    s += mv.getAdditions();
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    result += sign * mv.getResult() * getEntry(absoluteR, b); // adding sub-determinante times matrix entry
                                                                             // times appropriate sign
                    if (characteristic != 0) result = result % characteristic;
                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        s--; as--; // first addition was 0 + ..., so we do not count it
        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (iSB != 0) result = getReduction(result, iSB);
        IntMinorValue newMV(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
                                                           // number of retrievals does not make sense, as we
                                                           // do not use a cache.
        return newMV;
    }
}

/* This can only be used in the case of coefficients coming from a field!!! */
IntMinorValue IntMinorProcessor::getMinorPrivateBareiss(const int k, const MinorKey& mk, const int characteristic, const ideal& iSB) {
  assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
  int theRows[k]; mk.getAbsoluteRowIndices(theRows);
  int theColumns[k]; mk.getAbsoluteColumnIndices(theColumns);
  // the next line provides the return value for the case k = 1
  int e = getEntry(theRows[0], theColumns[0]);
  if (characteristic != 0) e = e % characteristic;
  if (iSB != 0) e = getReduction(e, iSB);
  IntMinorValue mv(e, 0, 0, 0, 0, -1, -1);
  if (k > 1)
  {
    // the matrix to perform Bareiss with
    long tempMatrix[k * k];
    // copy correct set of entries from _intMatrix to tempMatrix
    int i = 0;
    for (int r = 0; r < k; r++)
      for (int c = 0; c < k; c++)
      {
        e = getEntry(theRows[r], theColumns[c]);
        if (characteristic != 0) e = e % characteristic;
        tempMatrix[i++] = e;
      }
    // Bareiss algorithm operating on tempMatrix which is at least 2x2
    int sign = 1;   // This will store the correct sign resulting from permuting the rows of tempMatrix.
    int rowPermutation[k];   // This is for storing the permutation of rows resulting from searching for a
                             // non-zero pivot element.
    for (int i = 0; i < k; i++) rowPermutation[i] = i;
    int divisor = 1;   // the Bareiss divisor
    for (int r = 0; r <= k - 2; r++)
    {
      // look for a non-zero entry in column r:
      int i = r;
      while ((i < k) && (tempMatrix[rowPermutation[i] * k + r] == 0))
        i++;
      if (i == k)
        // There is no non-zero entry; hence the minor is zero.
        return IntMinorValue(0, 0, 0, 0, 0, -1, -1);
      if (i != r)
      {
        // We swap the rows with indices r and i:
        int j = rowPermutation[i];
        rowPermutation[i] = rowPermutation[r];
        rowPermutation[r] = j;
        // Now we know that tempMatrix[rowPermutation[r] * k + r] is not zero.
        // But carefull; we have to negate the sign, as there is always an odd
        // number of row transpositions to swap two given rows of a matrix.
        sign = -sign;
      }
      if (r >= 1) divisor = tempMatrix[rowPermutation[r - 1] * k + r - 1];
      for (int rr = r + 1; rr < k; rr++)
        for (int cc = r + 1; cc < k; cc++)
        {
          e = rowPermutation[rr] * k + cc;
          // Attention: The following may cause an overflow and thus a wrong result:
          tempMatrix[e] = tempMatrix[e] * tempMatrix[rowPermutation[r] * k + r]
                         - tempMatrix[rowPermutation[r] * k + cc] * tempMatrix[rowPermutation[rr] * k + r];
          // The following is, by theory, always a division without remainder:
          tempMatrix[e] = tempMatrix[e] / divisor;
          if (characteristic != 0) tempMatrix[e] = tempMatrix[e] % characteristic;
        }
    }
    int theValue = sign * tempMatrix[rowPermutation[k - 1] * k + k - 1];
    if (iSB != 0) theValue = getReduction(theValue, iSB);
    mv = IntMinorValue(theValue, 0, 0, 0, 0, -1, -1);
  }
  return mv;
}

int IntMinorProcessor::getEntry (const int rowIndex, const int columnIndex) const
{
  return _intMatrix[rowIndex * _columns + columnIndex];
}

IntMinorValue IntMinorProcessor::getMinorPrivateLaplace(const int k, const MinorKey& mk,
                                                        const bool multipleMinors,
                                                        Cache<MinorKey, IntMinorValue>& cch,
                                                        const int characteristic, const ideal& iSB) {
    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
    if (k == 1) {
        int e = getEntry(mk.getAbsoluteRowIndex(0), mk.getAbsoluteColumnIndex(0));
        if (characteristic != 0) e = e % characteristic;
        if (iSB != 0) e = getReduction(e, iSB);
        return IntMinorValue(e, 0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
    }
    else {
        int b = getBestLine(k, mk);                          // row or column with most zeros
        int result = 0;                                      // This will contain the value of the minor.
        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
                                                             // ..."a*" for accumulated operation counters
        IntMinorValue mv(0, 0, 0, 0, 0, 0, 0);               // for storing all intermediate minors
        if (b >= 0) {
            // This means that the best line is the row with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorRowKey:
            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
            for (int c = 0; c < k; c++) {
                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
                if (getEntry(b, absoluteC) != 0) { // Only then do we have to consider this sub-determinante.
                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is mk with row b and column absoluteC omitted.
                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
                        mv = cch.getValue(subMk);
                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
                                            // an impact on the internal ordering among cache entries.
                    }
                    else {
                        mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch, characteristic, iSB); // recursive call
                        // As this minor was not in the cache, we count the additions and
                        // multiplications that we needed to do in the recursive call:
                        m += mv.getMultiplications();
                        s += mv.getAdditions();
                    }
                    // In any case, we count all nested operations in the accumulative counters:
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    result += sign * mv.getResult() * getEntry(b, absoluteC); // adding sub-determinante times matrix entry
                                                                             // times appropriate sign
                    if (characteristic != 0) result = result % characteristic;
                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        else {
            b = - b - 1;
            // This means that the best line is the column with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorColumnKey:
            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
            for (int r = 0; r < k; r++) {
                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
                if (getEntry(absoluteR, b) != 0) { // Only then do we have to consider this sub-determinante.
                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is mk with row absoluteR and column b omitted.
                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
                        mv = cch.getValue(subMk);
                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
                                            // an impact on the internal ordering among cache entries.
                    }
                    else {
                        mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch, characteristic, iSB); // recursive call
                        // As this minor was not in the cache, we count the additions and
                        // multiplications that we needed to do in the recursive call:
                        m += mv.getMultiplications();
                        s += mv.getAdditions();
                    }
                    // In any case, we count all nested operations in the accumulative counters:
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    result += sign * mv.getResult() * getEntry(absoluteR, b); // adding sub-determinante times matrix entry
                                                                             // times appropriate sign
                    if (characteristic != 0) result = result % characteristic;
                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        // Let's cache the newly computed minor:
        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
        s--; as--; // first addition was 0 + ..., so we do not count it
        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (iSB != 0) result = getReduction(result, iSB);
        IntMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
        cch.put(mk, newMV); // Here's the actual put inside the cache.
        return newMV;
    }
}

PolyMinorProcessor::PolyMinorProcessor ()
{
  _polyMatrix = 0;
}

poly PolyMinorProcessor::getEntry (const int rowIndex, const int columnIndex) const
{
  return _polyMatrix[rowIndex * _columns + columnIndex];
}

bool PolyMinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
{
  return getEntry(absoluteRowIndex, absoluteColumnIndex) == NULL;
}

string PolyMinorProcessor::toString () const {
    char h[32];
    string t = "";
    string s = "PolyMinorProcessor:";
    s += "\n   matrix: ";
    sprintf(h, "%d", _rows); s += h;
    s += " x ";
    sprintf(h, "%d", _columns); s += h;
    int myIndexArray[500];
    s += "\n   considered submatrix has row indices: ";
    _container.getAbsoluteRowIndices(myIndexArray);
    for (int k = 0; k < _containerRows; k++) {
        if (k != 0) s += ", ";
        sprintf(h, "%d", myIndexArray[k]); s += h;
    }
    s += " (first row of matrix has index 0)";
    s += "\n   considered submatrix has column indices: ";
    _container.getAbsoluteColumnIndices(myIndexArray);
    for (int k = 0; k < _containerColumns; k++) {
        if (k != 0) s += ", ";
        sprintf(h, "%d", myIndexArray[k]); s += h;
    }
    s += " (first column of matrix has index 0)";
    s += "\n   size of considered minor(s): ";
    sprintf(h, "%d", _minorSize); s += h;
    s += "x";
    s += h;
    return s;
}

PolyMinorProcessor::~PolyMinorProcessor() {
    // free memory of _polyMatrix
    int n = _rows * _columns;
    for (int i = 0; i < n; i++)
      p_Delete(&_polyMatrix[i], currRing);
    delete [] _polyMatrix; _polyMatrix = 0;
}

void PolyMinorProcessor::defineMatrix (const int numberOfRows, const int numberOfColumns, const poly* polyMatrix) {
    // free memory of _polyMatrix
    int n = _rows * _columns;
    for (int i = 0; i < n; i++)
      p_Delete(&_polyMatrix[i], currRing);
    delete [] _polyMatrix; _polyMatrix = 0;

    _rows = numberOfRows;
    _columns = numberOfColumns;
    n = _rows * _columns;

    // allocate memory for new entries in _polyMatrix
    _polyMatrix = new poly[n];

    // copying values from one-dimensional method parameter "polyMatrix"
    for (int i = 0; i < n; i++)
      _polyMatrix[i] = pCopy(polyMatrix[i]);
}

PolyMinorValue PolyMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
                                            Cache<MinorKey, PolyMinorValue>& c, const ideal& iSB) {
    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
    _minorSize = dimension;
    // call a helper method which recursively computes the minor using the cache c:
    return getMinorPrivateLaplace(dimension, _container, false, c, iSB);
}

PolyMinorValue PolyMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
                                            const char* algorithm, const ideal& iSB) {
    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
    _minorSize = dimension;
    // call a helper method which computes the minor (without using a cache):
    if (strcmp(algorithm, "Laplace") == 0)
      return getMinorPrivateLaplace(_minorSize, _container, iSB);
    else if (strcmp(algorithm, "Bareiss") == 0)
      return getMinorPrivateBareiss(_minorSize, _container, iSB);
    else assume(false);
}

PolyMinorValue PolyMinorProcessor::getNextMinor(const char* algorithm, const ideal& iSB) {
    // call a helper method which computes the minor (without using a cache):
    if (strcmp(algorithm, "Laplace") == 0)
      return getMinorPrivateLaplace(_minorSize, _minor, iSB);
    else if (strcmp(algorithm, "Bareiss") == 0)
      return getMinorPrivateBareiss(_minorSize, _minor, iSB);
    else assume(false);
}

PolyMinorValue PolyMinorProcessor::getNextMinor(Cache<MinorKey, PolyMinorValue>& c, const ideal& iSB) {
    // computation with cache
    return getMinorPrivateLaplace(_minorSize, _minor, true, c, iSB);
}

/* assumes that all entries in polyMatrix are already reduced w.r.t. iSB */
PolyMinorValue PolyMinorProcessor::getMinorPrivateLaplace(const int k, const MinorKey& mk, const ideal& iSB) {
    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
    if (k == 1) {
      PolyMinorValue pmv(getEntry(mk.getAbsoluteRowIndex(0), mk.getAbsoluteColumnIndex(0)),
                         0, 0, 0, 0, -1, -1); // "-1" is to signal that any statistics about the
                                                 // number of retrievals does not make sense, as we do not use a cache.
      return pmv;
    }
    else {
        // Here, the minor must be 2x2 or larger.
        int b = getBestLine(k, mk);                          // row or column with most zeros
        poly result = NULL;                                  // This will contain the value of the minor.
        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
                                                             // ..."a*" for accumulated operation counters
        if (b >= 0) {
            // This means that the best line is the row with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorRowKey:
            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
            poly signPoly = NULL;
            for (int c = 0; c < k; c++) {
                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey with row b and column absoluteC omitted.
                    PolyMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, iSB);  // recursive call
                    m += mv.getMultiplications();
                    s += mv.getAdditions();
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    pDelete(&signPoly);
                    signPoly = pISet(sign);
                    poly temp = pp_Mult_qq(mv.getResult(), getEntry(b, absoluteC), currRing);
                    temp = p_Mult_q(signPoly, temp, currRing);
                    result = p_Add_q(result, temp, currRing);
                    signPoly = NULL;
                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        else {
            b = - b - 1;
            // This means that the best line is the column with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorColumnKey:
            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
            poly signPoly = NULL;
            for (int r = 0; r < k; r++) {
                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey with row absoluteR and column b omitted.
                    PolyMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, iSB);  // recursive call
                    m += mv.getMultiplications();
                    s += mv.getAdditions();
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    pDelete(&signPoly);
                    signPoly = pISet(sign);
                    poly temp = pp_Mult_qq(mv.getResult(), getEntry(absoluteR, b), currRing);
                    temp = p_Mult_q(signPoly, temp, currRing);
                    result = p_Add_q(result, temp, currRing);
                    signPoly = NULL;
                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        s--; as--; // first addition was 0 + ..., so we do not count it
        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (iSB != 0) result = kNF(iSB, currRing->qideal, result);
        PolyMinorValue newMV(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
                                                            // number of retrievals does not make sense, as we
                                                            // do not use a cache.
        pDelete(&result);
        return newMV;
    }
}

/* assumes that all entries in polyMatrix are already reduced w.r.t. iSB */
PolyMinorValue PolyMinorProcessor::getMinorPrivateLaplace(const int k, const MinorKey& mk,
                                                          const bool multipleMinors,
                                                          Cache<MinorKey, PolyMinorValue>& cch,
                                                          const ideal& iSB) {
    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
    if (k == 1) {
      PolyMinorValue pmv(getEntry(mk.getAbsoluteRowIndex(0), mk.getAbsoluteColumnIndex(0)),
                         0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
      return pmv;
    }
    else {
        int b = getBestLine(k, mk);                          // row or column with most zeros
        poly result = NULL;                                  // This will contain the value of the minor.
        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
                                                             // ..."a*" for accumulated operation counters
        if (b >= 0) {
            // This means that the best line is the row with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorRowKey:
            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
            poly signPoly = NULL;
            for (int c = 0; c < k; c++) {
                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
                    PolyMinorValue mv;              // for storing all intermediate minors
                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is mk with row b and column absoluteC omitted.
                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
                        mv = cch.getValue(subMk);
                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
                                            // an impact on the internal ordering among cache entries.
                    }
                    else {
                        mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch, iSB); // recursive call
                        // As this minor was not in the cache, we count the additions and
                        // multiplications that we needed to do in the recursive call:
                        m += mv.getMultiplications();
                        s += mv.getAdditions();
                    }
                    // In any case, we count all nested operations in the accumulative counters:
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    pDelete(&signPoly);
                    signPoly = pISet(sign);
                    poly temp = pp_Mult_qq(mv.getResult(), getEntry(b, absoluteC), currRing);
                    temp = p_Mult_q(signPoly, temp, currRing);
                    result = p_Add_q(result, temp, currRing);
                    signPoly = NULL;
                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        else {
            b = - b - 1;
            // This means that the best line is the column with absolute (0-based) index b.
            // Using Laplace, the sign of the contributing minors must be iterating;
            // the initial sign depends on the relative index of b in minorColumnKey:
            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
            poly signPoly = NULL;
            for (int r = 0; r < k; r++) {
                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
                    PolyMinorValue mv;              // for storing all intermediate minors
                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is mk with row absoluteR and column b omitted.
                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
                        mv = cch.getValue(subMk);
                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
                                            // an impact on the internal ordering among cache entries.
                    }
                    else {
                        mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch, iSB); // recursive call
                        // As this minor was not in the cache, we count the additions and
                        // multiplications that we needed to do in the recursive call:
                        m += mv.getMultiplications();
                        s += mv.getAdditions();
                    }
                    // In any case, we count all nested operations in the accumulative counters:
                    am += mv.getAccumulatedMultiplications();
                    as += mv.getAccumulatedAdditions();
                    pDelete(&signPoly);
                    signPoly = pISet(sign);
                    poly temp = pp_Mult_qq(mv.getResult(), getEntry(absoluteR, b), currRing);
                    temp = p_Mult_q(signPoly, temp, currRing);
                    result = p_Add_q(result, temp, currRing);
                    signPoly = NULL;
                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
                }
                sign = - sign; // alternating the sign
            }
        }
        // Let's cache the newly computed minor:
        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
        s--; as--; // first addition was 0 + ..., so we do not count it
        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
        if (iSB != 0) result = kNF(iSB, currRing->qideal, result);
        PolyMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
        pDelete(&result); result = NULL;
        cch.put(mk, newMV); // Here's the actual put inside the cache.
        return newMV;
    }
}

/* This can only be used in the case of coefficients coming from a field!!! */
void addOperationBucket(poly& f1, poly& f2, kBucket_pt& bucket)
{
  /* fills all terms of f1 * f2 into the bucket */
  poly a = f1; poly b = f2;
  int aLen = pLength(a); int bLen = pLength(b);
  if (aLen > bLen)
  {
    b = f1; a = f2; bLen = aLen;
  }
  pNormalize(b);

  while (a != NULL)
  {
    /* The next line actually uses only LT(a): */
    kBucket_Plus_mm_Mult_pp(bucket, a, b, bLen);
    a = pNext(a);
  }
}

/* computes the polynomial (p1 * p2 - p3 * p4) and puts result into p1;
   the method destroys the old value of p1;
   p2, p3, and p4 may be pNormalize-d but must, apart from that,
   not be changed;
   This can only be used in the case of coefficients coming from a field!!! */
void elimOperationBucketNoDiv(poly &p1, poly &p2, poly &p3, poly &p4)
{
  kBucket_pt myBucket = kBucketCreate();
  addOperationBucket(p1, p2, myBucket);
  poly p3Neg = pNeg(pCopy(p3));
  addOperationBucket(p3Neg, p4, myBucket);
  pDelete(&p3Neg);
  pDelete(&p1);
  p1 = kBucketClear(myBucket);
  kBucketDestroy(&myBucket);
}

/* computes the polynomial (p1 * p2 - p3 * p4) / p5 and puts result into p1;
   the method destroys the old value of p1;
   p2, p3, p4, and p5 may be pNormalize-d but must, apart from that,
   not be changed;
   c5 is assumed to be the leading coefficient of p5;
   p5Len is assumed to be the length of p5;
   This can only be used in the case of coefficients coming from a field!!! */
void elimOperationBucket(poly &p1, poly &p2, poly &p3, poly &p4, poly &p5, number &c5, int p5Len)
{
  kBucket_pt myBucket = kBucketCreate();
  addOperationBucket(p1, p2, myBucket);
  poly p3Neg = pNeg(pCopy(p3));
  addOperationBucket(p3Neg, p4, myBucket);
  pDelete(&p3Neg);

  // Now, myBucket contains all terms of p1 * p2 - p3 * p4.
  // Now we need to perform the polynomial division myBucket / p5
  // which is known to work without remainder:
  pDelete(&p1); poly helperPoly = NULL;

  poly bucketLm = pCopy(kBucketGetLm(myBucket));
  while (bucketLm != NULL)
  {
    /* divide bucketLm by the leading term of p5 and put result into bucketLm;
       we start with the coefficients;
       note that bucketLm will always represent a term */
    number coeff = nDiv(pGetCoeff(bucketLm), c5);
    nNormalize(coeff);
    pSetCoeff(bucketLm, coeff);
    /* subtract exponent vector of p5 from that of quotient; modifies quotient */
    p_ExpVectorSub(bucketLm, p5, currRing);
    kBucket_Minus_m_Mult_p(myBucket, bucketLm, p5, &p5Len);
    /* The following lines make bucketLm the new leading term of p1,
       i.e., put bucketLm in front of everything which is already in p1.
       Thus, after the while loop, we need to revert p1. */
    helperPoly = bucketLm;
    helperPoly->next = p1;
    p1 = helperPoly;

    bucketLm = pCopy(kBucketGetLm(myBucket));
  }
  p1 = pReverse(p1);
  kBucketDestroy(&myBucket);
}

/* assumes that all entries in polyMatrix are already reduced w.r.t. iSB
   This can only be used in the case of coefficients coming from a field!!! */
PolyMinorValue PolyMinorProcessor::getMinorPrivateBareiss(const int k, const MinorKey& mk, const ideal& iSB) {
  assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
  int theRows[k]; mk.getAbsoluteRowIndices(theRows);
  int theColumns[k]; mk.getAbsoluteColumnIndices(theColumns);
  if (k == 1)
    return PolyMinorValue(getEntry(theRows[0], theColumns[0]), 0, 0, 0, 0, -1, -1);
  else /* k > 0 */
  {
    // the matrix to perform Bareiss with
    poly* tempMatrix = (poly*)omAlloc(k * k * sizeof(poly));
    // copy correct set of entries from _polyMatrix to tempMatrix
    int i = 0;
    for (int r = 0; r < k; r++)
      for (int c = 0; c < k; c++)
        tempMatrix[i++] = pCopy(getEntry(theRows[r], theColumns[c]));

    // Bareiss algorithm operating on tempMatrix which is at least 2x2
    int sign = 1;   // This will store the correct sign resulting from permuting the rows of tempMatrix.
    int rowPermutation[k];   // This is for storing the permutation of rows resulting from searching for a
                             // non-zero pivot element.
    for (int i = 0; i < k; i++) rowPermutation[i] = i;
    poly divisor = NULL;
    int divisorLength = 0;
    number divisorLC;
    for (int r = 0; r <= k - 2; r++)
    {
      /* look for a non-zero entry in column r, rows = r .. (k - 1)
         s.t. the polynomial has least complexity: */
      int minComplexity = -1; int complexity = 0; int bestRow = -1; poly pp = NULL;
      for (int i = r; i < k; i++)
      {
        pp = tempMatrix[rowPermutation[i] * k + r];
        if (pp != NULL)
        {
          if (minComplexity == -1)
          {
            minComplexity = pSize(pp);
            bestRow = i;
          }
          else
          {
            complexity = 0;
            while ((pp != NULL) && (complexity < minComplexity)) { complexity += nSize(pGetCoeff(pp)); pp = pNext(pp); }
            if (complexity < minComplexity)
            {
              minComplexity = complexity;
              bestRow = i;
            }
          }
          if (minComplexity <= 1) break; /* terminate the search */
        }
      }
      if (bestRow == -1)
      {
        // There is no non-zero entry; hence the minor is zero.
        for (int i = 0; i < k * k; i++) pDelete(&tempMatrix[i]);
        return PolyMinorValue(NULL, 0, 0, 0, 0, -1, -1);
      }
      pNormalize(tempMatrix[rowPermutation[bestRow] * k + r]);
      if (bestRow != r)
      {
        // We swap the rows with indices r and i:
        int j = rowPermutation[bestRow];
        rowPermutation[bestRow] = rowPermutation[r];
        rowPermutation[r] = j;
        // Now we know that tempMatrix[rowPermutation[r] * k + r] is not zero.
        // But carefull; we have to negate the sign, as there is always an odd
        // number of row transpositions to swap two given rows of a matrix.
        sign = -sign;
      }
      if (r != 0)
      {
        divisor = tempMatrix[rowPermutation[r - 1] * k + r - 1];
        pNormalize(divisor);
        divisorLength = pLength(divisor);
        divisorLC = pGetCoeff(divisor);
      }
      for (int rr = r + 1; rr < k; rr++)
        for (int cc = r + 1; cc < k; cc++)
        {
          if (r == 0)
            elimOperationBucketNoDiv(tempMatrix[rowPermutation[rr] * k + cc],
                                     tempMatrix[rowPermutation[r]  * k + r],
                                     tempMatrix[rowPermutation[r]  * k + cc],
                                     tempMatrix[rowPermutation[rr] * k + r]);
          else
            elimOperationBucket(tempMatrix[rowPermutation[rr] * k + cc],
                                tempMatrix[rowPermutation[r]  * k + r],
                                tempMatrix[rowPermutation[r]  * k + cc],
                                tempMatrix[rowPermutation[rr] * k + r],
                                divisor, divisorLC, divisorLength);
        }
    }
    poly result = tempMatrix[rowPermutation[k - 1] * k + k - 1];
    if (sign == -1) result = pNeg(result);
    if (iSB != 0) result = kNF(iSB, currRing->qideal, result);
    PolyMinorValue mv(result, 0, 0, 0, 0, -1, -1);
    for (int i = 0; i < k * k; i++) pDelete(&tempMatrix[i]);
    omFreeSize(tempMatrix, k * k * sizeof(poly));
    return mv;
  }
}