source: git/Singular/MinorProcessor.cc @ 83f349c

spielwiese
Last change on this file since 83f349c was 83f349c, checked in by Frank Seelisch <seelisch@…>, 14 years ago
git-svn-id: file:///usr/local/Singular/svn/trunk@12250 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 40.3 KB
Line 
1#include "mod2.h"
2
3#ifdef HAVE_MINOR
4
5#include "structs.h"
6#include "polys.h"
7#include <MinorProcessor.h>
8#include "febase.h"
9
10void MinorProcessor::print() const {
11    PrintS(this->toString().c_str());
12}
13
14int MinorProcessor::getBestLine (const int k, const MinorKey& mk) const {
15    // This method identifies the row or column with the most zeros.
16    // The returned index (bestIndex) is absolute within the pre-defined matrix.
17    // If some row has the most zeros, then the absolute (0-based) row index is returned.
18    // If, contrariwise, some column has the most zeros, then -1 minus the absolute (0-based) column index is returned.
19    int numberOfZeros = 0;
20    int bestIndex = 100000;    // We start with an invalid row/column index.
21    int maxNumberOfZeros = -1; // We update this variable whenever we find a new so-far optimal row or column.
22    for (int r = 0; r < k; r++) {
23        // iterate through all k rows of the momentary minor
24        int absoluteR = mk.getAbsoluteRowIndex(r);
25        numberOfZeros = 0;
26        for (int c = 0; c < k; c++) {
27            int absoluteC = mk.getAbsoluteColumnIndex(c);
28            if (isEntryZero(absoluteR, absoluteC)) numberOfZeros++;
29        }
30        if (numberOfZeros > maxNumberOfZeros) {
31            // We found a new best line which is a row.
32            bestIndex = absoluteR;
33            maxNumberOfZeros = numberOfZeros;
34        }
35    };
36    for (int c = 0; c < k; c++) {
37        int absoluteC = mk.getAbsoluteColumnIndex(c);
38        numberOfZeros = 0;
39        for (int r = 0; r < k; r++) {
40            int absoluteR = mk.getAbsoluteRowIndex(r);
41            if (isEntryZero(absoluteR, absoluteC)) numberOfZeros++;
42        }
43        if (numberOfZeros > maxNumberOfZeros) {
44            // We found a new best line which is a column. So we transform the return value.
45            // Note that we can easily get back absoluteC from bestLine: absoluteC = - 1 - bestLine.
46            bestIndex = - absoluteC - 1;
47            maxNumberOfZeros = numberOfZeros;
48        }
49    };
50    return bestIndex;
51}
52
53void MinorProcessor::setMinorSize(const int minorSize) {
54    _minorSize = minorSize;
55    _minor.reset();
56}
57
58bool MinorProcessor::hasNextMinor() {
59    return setNextKeys(_minorSize);
60}
61
62void MinorProcessor::getCurrentRowIndices(int* const target) const {
63    return _minor.getAbsoluteRowIndices(target);
64}
65
66void MinorProcessor::getCurrentColumnIndices(int* const target) const {
67    return _minor.getAbsoluteColumnIndices(target);
68}
69
70void MinorProcessor::defineSubMatrix(const int numberOfRows, const int* rowIndices,
71                                     const int numberOfColumns, const int* columnIndices) {
72    // The method assumes ascending row and column indices in the two argument arrays.
73    // These indices are understood to be zero-based.
74    // The method will set the two arrays of ints in _container.
75    // Example: The indices 0, 2, 3, 7 will be converted to an array with one int
76    // representing the binary number 10001101 (check bits from right to left).
77
78    _containerRows = numberOfRows;
79    int highestRowIndex = rowIndices[numberOfRows - 1];
80    int rowBlockCount = (highestRowIndex / 32) + 1;
81    unsigned int rowBlocks[rowBlockCount];
82    for (int i = 0; i < rowBlockCount; i++) rowBlocks[i] = 0;
83    for (int i = 0; i < numberOfRows; i++) {
84        int blockIndex = rowIndices[i] / 32;
85        int offset = rowIndices[i] % 32;
86        rowBlocks[blockIndex] += (1 << offset);
87    }
88
89    _containerColumns = numberOfColumns;
90    int highestColumnIndex = columnIndices[numberOfColumns - 1];
91    int columnBlockCount = (highestColumnIndex / 32) + 1;
92    unsigned int columnBlocks[columnBlockCount];
93    for (int i = 0; i < columnBlockCount; i++) columnBlocks[i] = 0;
94    for (int i = 0; i < numberOfColumns; i++) {
95        int blockIndex = columnIndices[i] / 32;
96        int offset = columnIndices[i] % 32;
97        columnBlocks[blockIndex] += (1 << offset);
98    }
99
100    _container.set(rowBlockCount, rowBlocks, columnBlockCount, columnBlocks);
101}
102
103bool MinorProcessor::setNextKeys(const int k) {
104    // This method moves _minor to the next valid kxk-minor within _container.
105    // It returns true iff this is successful, i.e. iff _minor did not already encode the final kxk-minor.
106    if (_minor.compare(MinorKey(0, 0, 0, 0)) == 0) {
107        // This means that we haven't started yet. Thus, we are about to compute the first kxk-minor.
108        _minor.selectFirstRows(k, _container);
109        _minor.selectFirstColumns(k, _container);
110        return true;
111    }
112    else if (_minor.selectNextColumns(k, _container)) {
113        // Here we were able to pick a next subset of columns (within the same subset of rows).
114        return true;
115    }
116    else if (_minor.selectNextRows(k, _container)) {
117        // Here we were not able to pick a next subset of columns (within the same subset of rows).
118        // But we could pick a next subset of rows. We must hence reset the subset of columns:
119        _minor.selectFirstColumns(k, _container);
120        return true;
121    }
122    else {
123        // We were neither able to pick a next subset of columns nor of rows.
124        // I.e., we have iterated through all sensible choices of subsets of rows and columns.
125        return false;
126    }
127}
128
129bool MinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
130{
131  assume(false);
132  return false;
133}
134
135string MinorProcessor::toString () const
136{
137  assume(false);
138  return "";
139}
140
141int MinorProcessor::IOverJ(const int i, const int j) {
142    // This is a non-recursive implementation.
143    assert(i >= 0 && j >= 0 && i >= j);
144    if (j == 0 || i == j) return 1;
145    int result = 1;
146    for (int k = i - j + 1; k <= i; k++) result *= k;
147    // Here, we have result = (i - j + 1) * ... * i.
148    for (int k = 2; k <= j; k++) result /= k;
149    // Here, we have result = (i - j + 1) * ... * i / 1 / 2 ... / j = i! / j! / (i - j)!.
150    return result;
151}
152
153int MinorProcessor::Faculty(const int i) {
154    // This is a non-recursive implementation.
155    assert(i >= 0);
156    int result = 1;
157    for (int j = 1; j <= i; j++) result *= j;
158    // Here, we have result = 1 * 2 * ... * i = i!
159    return result;
160}
161
162int MinorProcessor::NumberOfRetrievals (const int rows, const int columns, const int containerMinorSize,
163                        const int minorSize, const bool multipleMinors) {
164    // This method computes the number of potential retrievals of a single minor when computing
165    // all minors of a given size within a matrix of given size.
166    int result = 0;
167    if (multipleMinors) {
168        // Here, we would like to compute all minors of size
169        // containerMinorSize x containerMinorSize in a matrix of size rows x columns.
170        // Then, we need to retrieve any minor of size minorSize x minorSize exactly
171        // n times, where n is as follows:
172        result = IOverJ(rows - minorSize, containerMinorSize - minorSize)
173               * IOverJ(columns - minorSize, containerMinorSize - minorSize)
174               * Faculty(containerMinorSize - minorSize);
175    }
176    else {
177        // Here, we would like to compute just one minor of size
178        // containerMinorSize x containerMinorSize. Then, we need to retrieve
179        // any minor of size minorSize x minorSize exactly
180        // (containerMinorSize - minorSize)! times:
181        result = Faculty(containerMinorSize - minorSize);
182    }
183    return result;
184}
185
186MinorProcessor::MinorProcessor () {
187    _container = MinorKey(0, 0, 0, 0);
188    _minor = MinorKey(0, 0, 0, 0);
189    _containerRows = 0;
190    _containerColumns = 0;
191    _minorSize = 0;
192    _rows = 0;
193    _columns = 0;
194}
195
196IntMinorProcessor::IntMinorProcessor () {
197    _matrix = 0;
198}
199
200string IntMinorProcessor::toString () const {
201    char h[32];
202    string t = "";
203    string s = "IntMinorProcessor:";
204    s += "\n   matrix: ";
205    sprintf(h, "%d", _rows); s += h;
206    s += " x ";
207    sprintf(h, "%d", _columns); s += h;
208    for (int r = 0; r < _rows; r++) {
209        s += "\n      ";
210        for (int c = 0; c < _columns; c++) {
211            sprintf(h, "%d", _matrix[r][c]); t = h;
212            for (int k = 0; k < int(4 - strlen(h)); k++) s += " ";
213            s += t;
214        }
215    }
216    int myIndexArray[500];
217    s += "\n   considered submatrix has row indices: ";
218    _container.getAbsoluteRowIndices(myIndexArray);
219    for (int k = 0; k < _containerRows; k++) {
220        if (k != 0) s += ", ";
221        sprintf(h, "%d", myIndexArray[k]); s += h;
222    }
223    s += " (first row of matrix has index 0)";
224    s += "\n   considered submatrix has column indices: ";
225    _container.getAbsoluteColumnIndices(myIndexArray);
226    for (int k = 0; k < _containerColumns; k++) {
227        if (k != 0) s += ", ";
228        sprintf(h, "%d", myIndexArray[k]); s += h;
229    }
230    s += " (first column of matrix has index 0)";
231    s += "\n   size of considered minor(s): ";
232    sprintf(h, "%d", _minorSize); s += h;
233    s += "x";
234    s += h;
235    return s;
236}
237
238bool IntMinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
239{
240  return _matrix[absoluteRowIndex][absoluteColumnIndex] == 0;
241}
242
243IntMinorProcessor::~IntMinorProcessor() {
244    // free memory of _matrix
245    for (int i = 0; i < _rows; i++) {
246        delete [] _matrix[i]; _matrix[i] = 0;
247    }
248    delete [] _matrix; _matrix = 0;
249}
250
251void IntMinorProcessor::defineMatrix (const int numberOfRows, const int numberOfColumns, const int* matrix) {
252    // free memory of _matrix
253    for (int i = 0; i < _rows; i++) {
254        delete [] _matrix[i]; _matrix[i] = 0;
255    }
256    delete [] _matrix; _matrix = 0;
257
258    _rows = numberOfRows;
259    _columns = numberOfColumns;
260
261    // allocate memory for new entries in _matrix
262    _matrix = new int*[_rows];
263    for (int i = 0; i < _rows; i++) _matrix[i] = new int[_columns];
264
265    // copying values from one-dimensional method parameter "matrix"
266    for (int r = 0; r < _rows; r++)
267        for (int c = 0; c < _columns; c++)
268            _matrix[r][c] = matrix[r * _columns + c];
269}
270
271IntMinorValue IntMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
272                                          Cache<MinorKey, IntMinorValue>& c, const int characteristic) {
273    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
274    _minorSize = dimension;
275    // call a helper method which recursively computes the minor using the cache c:
276    return getMinorPrivate(dimension, _container, false, c, characteristic);
277}
278
279IntMinorValue IntMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices, const int characteristic) {
280    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
281    _minorSize = dimension;
282    // call a helper method which recursively computes the minor (without using a cache):
283    return getMinorPrivate(_minorSize, _container, characteristic);
284}
285
286IntMinorValue IntMinorProcessor::getNextMinor(const int characteristic) {
287    // computation without cache
288    return getMinorPrivate(_minorSize, _minor, characteristic);
289}
290
291IntMinorValue IntMinorProcessor::getNextMinor(Cache<MinorKey, IntMinorValue>& c, const int characteristic) {
292    // computation with cache
293    return getMinorPrivate(_minorSize, _minor, true, c, characteristic);
294}
295
296IntMinorValue IntMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk, const int characteristic) {
297    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
298    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
299    if (k == 1) {
300        return IntMinorValue(_matrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
301                              0, 0, 0, 0, -1, -1); // "-1" is to signal that any statistics about the
302                                                   // number of retrievals does not make sense, as we do not use a cache.
303    }
304    else {
305        // Here, the minor must be 2x2 or larger.
306        int b = getBestLine(k, mk);                          // row or column with most zeros
307        int result = 0;                                      // This will contain the value of the minor.
308        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
309                                                             // ..."a*" for accumulated operation counters
310        if (b >= 0) {
311            // This means that the best line is the row with absolute (0-based) index b.
312            // Using Laplace, the sign of the contributing minors must be iterating;
313            // the initial sign depends on the relative index of b in minorRowKey:
314            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
315            for (int c = 0; c < k; c++) {
316                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
317                if (_matrix[b][absoluteC] != 0) { // Only then do we have to consider this sub-determinante.
318                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is mk with row b and column absoluteC omitted.
319                    IntMinorValue mv = getMinorPrivate(k - 1, subMk, characteristic);  // recursive call
320                    m += mv.getMultiplications();
321                    s += mv.getAdditions();
322                    am += mv.getAccumulatedMultiplications();
323                    as += mv.getAccumulatedAdditions();
324                    result += sign * mv.getResult() * _matrix[b][absoluteC]; // adding sub-determinante times matrix entry
325                                                                             // times appropriate sign
326                    if (characteristic != 0) result = result % characteristic;
327                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
328                }
329                sign = - sign; // alternating the sign
330            }
331        }
332        else {
333            b = - b - 1;
334            // This means that the best line is the column with absolute (0-based) index b.
335            // Using Laplace, the sign of the contributing minors must be iterating;
336            // the initial sign depends on the relative index of b in minorColumnKey:
337            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
338            for (int r = 0; r < k; r++) {
339                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
340                if (_matrix[absoluteR][b] != 0) { // Only then do we have to consider this sub-determinante.
341                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is mk with row absoluteR and column b omitted.
342                    IntMinorValue mv = getMinorPrivate(k - 1, subMk, characteristic);  // recursive call
343                    m += mv.getMultiplications();
344                    s += mv.getAdditions();
345                    am += mv.getAccumulatedMultiplications();
346                    as += mv.getAccumulatedAdditions();
347                    result += sign * mv.getResult() * _matrix[absoluteR][b]; // adding sub-determinante times matrix entry
348                                                                             // times appropriate sign
349                    if (characteristic != 0) result = result % characteristic;
350                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
351                }
352                sign = - sign; // alternating the sign
353            }
354        }
355        s--; as--; // first addition was 0 + ..., so we do not count it
356        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
357        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
358        IntMinorValue newMV(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
359                                                           // number of retrievals does not make sense, as we
360                                                           // do not use a cache.
361        return newMV;
362    }
363}
364
365IntMinorValue IntMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk,
366                                                 const bool multipleMinors,
367                                                 Cache<MinorKey, IntMinorValue>& cch,
368                                                 const int characteristic) {
369    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
370    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
371    if (k == 1) {
372        return IntMinorValue(_matrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
373                              0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
374    }
375    else {
376        int b = getBestLine(k, mk);                          // row or column with most zeros
377        int result = 0;                                      // This will contain the value of the minor.
378        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
379                                                             // ..."a*" for accumulated operation counters
380        IntMinorValue mv(0, 0, 0, 0, 0, 0, 0);               // for storing all intermediate minors
381        if (b >= 0) {
382            // This means that the best line is the row with absolute (0-based) index b.
383            // Using Laplace, the sign of the contributing minors must be iterating;
384            // the initial sign depends on the relative index of b in minorRowKey:
385            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
386            for (int c = 0; c < k; c++) {
387                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
388                if (_matrix[b][absoluteC] != 0) { // Only then do we have to consider this sub-determinante.
389                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is mk with row b and column absoluteC omitted.
390                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
391                        mv = cch.getValue(subMk);
392                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
393                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
394                                            // an impact on the internal ordering among cache entries.
395                    }
396                    else {
397                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch, characteristic); // recursive call
398                        // As this minor was not in the cache, we count the additions and
399                        // multiplications that we needed to do in the recursive call:
400                        m += mv.getMultiplications();
401                        s += mv.getAdditions();
402                    }
403                    // In any case, we count all nested operations in the accumulative counters:
404                    am += mv.getAccumulatedMultiplications();
405                    as += mv.getAccumulatedAdditions();
406                    result += sign * mv.getResult() * _matrix[b][absoluteC]; // adding sub-determinante times matrix entry
407                                                                             // times appropriate sign
408                    if (characteristic != 0) result = result % characteristic;
409                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
410                }
411                sign = - sign; // alternating the sign
412            }
413        }
414        else {
415            b = - b - 1;
416            // This means that the best line is the column with absolute (0-based) index b.
417            // Using Laplace, the sign of the contributing minors must be iterating;
418            // the initial sign depends on the relative index of b in minorColumnKey:
419            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
420            for (int r = 0; r < k; r++) {
421                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
422                if (_matrix[absoluteR][b] != 0) { // Only then do we have to consider this sub-determinante.
423                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is mk with row absoluteR and column b omitted.
424                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
425                        mv = cch.getValue(subMk);
426                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
427                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
428                                            // an impact on the internal ordering among cache entries.
429                    }
430                    else {
431                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch, characteristic); // recursive call
432                        // As this minor was not in the cache, we count the additions and
433                        // multiplications that we needed to do in the recursive call:
434                        m += mv.getMultiplications();
435                        s += mv.getAdditions();
436                    }
437                    // In any case, we count all nested operations in the accumulative counters:
438                    am += mv.getAccumulatedMultiplications();
439                    as += mv.getAccumulatedAdditions();
440                    result += sign * mv.getResult() * _matrix[absoluteR][b]; // adding sub-determinante times matrix entry
441                                                                             // times appropriate sign
442                    if (characteristic != 0) result = result % characteristic;
443                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
444                }
445                sign = - sign; // alternating the sign
446            }
447        }
448        // Let's cache the newly computed minor:
449        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
450        s--; as--; // first addition was 0 + ..., so we do not count it
451        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
452        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
453        IntMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
454        cch.put(mk, newMV); // Here's the actual put inside the cache.
455        return newMV;
456    }
457}
458
459PolyMinorProcessor::PolyMinorProcessor () {
460    _polyMatrix = 0;
461}
462
463string PolyMinorProcessor::toString () const {
464    char h[32];
465    string t = "";
466    string s = "PolyMinorProcessor:";
467    s += "\n   matrix: ";
468    sprintf(h, "%d", _rows); s += h;
469    s += " x ";
470    sprintf(h, "%d", _columns); s += h;
471    int myIndexArray[500];
472    s += "\n   considered submatrix has row indices: ";
473    _container.getAbsoluteRowIndices(myIndexArray);
474    for (int k = 0; k < _containerRows; k++) {
475        if (k != 0) s += ", ";
476        sprintf(h, "%d", myIndexArray[k]); s += h;
477    }
478    s += " (first row of matrix has index 0)";
479    s += "\n   considered submatrix has column indices: ";
480    _container.getAbsoluteColumnIndices(myIndexArray);
481    for (int k = 0; k < _containerColumns; k++) {
482        if (k != 0) s += ", ";
483        sprintf(h, "%d", myIndexArray[k]); s += h;
484    }
485    s += " (first column of matrix has index 0)";
486    s += "\n   size of considered minor(s): ";
487    sprintf(h, "%d", _minorSize); s += h;
488    s += "x";
489    s += h;
490    return s;
491}
492
493bool PolyMinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
494{
495  return (_polyMatrix[absoluteRowIndex][absoluteColumnIndex] == NULL);
496}
497
498PolyMinorProcessor::~PolyMinorProcessor() {
499    // free memory of _polyMatrix
500    for (int i = 0; i < _rows; i++) {
501        for (int j = 0; j < _columns; j++) {
502            p_Delete(&_polyMatrix[i][j], currRing);
503        }
504        delete [] _polyMatrix[i]; _polyMatrix[i] = 0;
505    }
506    delete [] _polyMatrix; _polyMatrix = 0;
507}
508
509void PolyMinorProcessor::defineMatrix (const int numberOfRows, const int numberOfColumns, const poly* polyMatrix) {
510    // free memory of _polyMatrix
511    for (int i = 0; i < _rows; i++) {
512        for (int j = 0; j < _columns; j++)
513            p_Delete(&_polyMatrix[i][j], currRing);
514        delete [] _polyMatrix[i]; _polyMatrix[i] = 0;
515    }
516    delete [] _polyMatrix; _polyMatrix = 0;
517
518    _rows = numberOfRows;
519    _columns = numberOfColumns;
520
521    // allocate memory for new entries in _matrix
522    _polyMatrix = new poly*[_rows];
523    for (int i = 0; i < _rows; i++) _polyMatrix[i] = new poly[_columns];
524
525    // copying values from one-dimensional method parameter "matrix"
526    for (int r = 0; r < _rows; r++)
527        for (int c = 0; c < _columns; c++)
528            _polyMatrix[r][c] = pCopy(polyMatrix[r * _columns + c]);
529}
530
531PolyMinorValue PolyMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
532                                            Cache<MinorKey, PolyMinorValue>& c) {
533    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
534    _minorSize = dimension;
535    // call a helper method which recursively computes the minor using the cache c:
536    return getMinorPrivate(dimension, _container, false, c);
537}
538
539PolyMinorValue PolyMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices) {
540    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
541    _minorSize = dimension;
542    // call a helper method which recursively computes the minor (without using a cache):
543    return getMinorPrivate(_minorSize, _container);
544}
545
546PolyMinorValue PolyMinorProcessor::getNextMinor() {
547    // computation without cache
548    return getMinorPrivate(_minorSize, _minor);
549}
550
551PolyMinorValue PolyMinorProcessor::getNextMinor(Cache<MinorKey, PolyMinorValue>& c) {
552    // computation with cache
553    return getMinorPrivate(_minorSize, _minor, true, c);
554}
555
556PolyMinorValue PolyMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk) {
557    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
558    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
559    if (k == 1) {
560        return PolyMinorValue(_polyMatrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
561                              0, 0, 0, 0, -1, -1); // "-1" is to signal that any statistics about the
562                                                   // number of retrievals does not make sense, as we do not use a cache.
563    }
564    else {
565        // Here, the minor must be 2x2 or larger.
566        int b = getBestLine(k, mk);                          // row or column with most zeros
567        poly result = NULL;                                  // This will contain the value of the minor.
568        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
569                                                             // ..."a*" for accumulated operation counters
570        if (b >= 0) {
571            // This means that the best line is the row with absolute (0-based) index b.
572            // Using Laplace, the sign of the contributing minors must be iterating;
573            // the initial sign depends on the relative index of b in minorRowKey:
574            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
575            poly signPoly = NULL;
576            for (int c = 0; c < k; c++) {
577                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
578                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
579                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey with row b and column absoluteC omitted.
580                    PolyMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
581                    m += mv.getMultiplications();
582                    s += mv.getAdditions();
583                    am += mv.getAccumulatedMultiplications();
584                    as += mv.getAccumulatedAdditions();
585                    pDelete(&signPoly);
586                    signPoly = pISet(sign);
587                    poly temp = pp_Mult_qq(mv.getResult(), _polyMatrix[b][absoluteC], currRing);
588                    temp = p_Mult_q(signPoly, temp, currRing);
589                    result = p_Add_q(result, temp, currRing);
590                    signPoly = NULL;
591                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
592                }
593                sign = - sign; // alternating the sign
594            }
595        }
596        else {
597            b = - b - 1;
598            // This means that the best line is the column with absolute (0-based) index b.
599            // Using Laplace, the sign of the contributing minors must be iterating;
600            // the initial sign depends on the relative index of b in minorColumnKey:
601            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
602            poly signPoly = NULL;
603            for (int r = 0; r < k; r++) {
604                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
605                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
606                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey with row absoluteR and column b omitted.
607                    PolyMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
608                    m += mv.getMultiplications();
609                    s += mv.getAdditions();
610                    am += mv.getAccumulatedMultiplications();
611                    as += mv.getAccumulatedAdditions();
612                    pDelete(&signPoly);
613                    signPoly = pISet(sign);
614                    poly temp = pp_Mult_qq(mv.getResult(), _polyMatrix[absoluteR][b], currRing);
615                    temp = p_Mult_q(signPoly, temp, currRing);
616                    result = p_Add_q(result, temp, currRing);
617                    signPoly = NULL;
618                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
619                }
620                sign = - sign; // alternating the sign
621            }
622        }
623        s--; as--; // first addition was 0 + ..., so we do not count it
624        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
625        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
626        PolyMinorValue newMV(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
627                                                            // number of retrievals does not make sense, as we
628                                                            // do not use a cache.
629//printf("\nMINOR to put: %s", pString(result));
630        pDelete(&result);
631        return newMV;
632    }
633}
634
635PolyMinorValue PolyMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk,
636                                                   const bool multipleMinors,
637                                                   Cache<MinorKey, PolyMinorValue>& cch) {
638//omUpdateInfo(); printf("\nused bytes: %12ld, called: %s  (k = %d)", om_Info.UsedBytes, "getMinorPrivate", k);
639    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
640    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
641    if (k == 1) {
642        return PolyMinorValue(_polyMatrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
643                              0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
644    }
645    else {
646        int b = getBestLine(k, mk);                          // row or column with most zeros
647        poly result = NULL;                                  // This will contain the value of the minor.
648//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d)", om_Info.UsedBytes, "poly result = NULL;", k);
649        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
650                                                             // ..."a*" for accumulated operation counters
651        if (b >= 0) {
652            // This means that the best line is the row with absolute (0-based) index b.
653            // Using Laplace, the sign of the contributing minors must be iterating;
654            // the initial sign depends on the relative index of b in minorRowKey:
655            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
656            poly signPoly = NULL;
657            for (int c = 0; c < k; c++) {
658                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
659//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "MinorKey subMk = mk.getSubMinorKey(b, absoluteC);", k);
660                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
661                    PolyMinorValue mv;              // for storing all intermediate minors
662                    MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is mk with row b and column absoluteC omitted.
663//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "PolyMinorValue mv;", k);
664                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
665                        mv = cch.getValue(subMk);
666                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
667                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
668                                            // an impact on the internal ordering among cache entries.
669//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "cch.put(subMk, mv);", k);
670                    }
671                    else {
672//omUpdateInfo(); printf("\nused bytes: %12ld, before: %s  (k = %d, case 1)", om_Info.UsedBytes, "mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch);", k);
673                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
674//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch);", k);
675                        // As this minor was not in the cache, we count the additions and
676                        // multiplications that we needed to do in the recursive call:
677                        m += mv.getMultiplications();
678                        s += mv.getAdditions();
679                    }
680                    // In any case, we count all nested operations in the accumulative counters:
681                    am += mv.getAccumulatedMultiplications();
682                    as += mv.getAccumulatedAdditions();
683                    pDelete(&signPoly);
684//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "pDelete(&signPoly);", k);
685                    signPoly = pISet(sign);
686//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "signPoly = pISet(sign);", k);
687//printf("\nmv.getResult() = %s", pString(mv.getResult()));
688                    poly temp = pp_Mult_qq(mv.getResult(), _polyMatrix[b][absoluteC], currRing);
689                    temp = p_Mult_q(signPoly, temp, currRing);
690                    result = p_Add_q(result, temp, currRing);
691//printf("\nresult = %s", pString(result));
692//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 1)", om_Info.UsedBytes, "ops(result, signPoly, mv.getResult(), _polyMatrix[b][absoluteC]);", k);
693                    signPoly = NULL;
694                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
695                }
696                sign = - sign; // alternating the sign
697            }
698        }
699        else {
700            b = - b - 1;
701            // This means that the best line is the column with absolute (0-based) index b.
702            // Using Laplace, the sign of the contributing minors must be iterating;
703            // the initial sign depends on the relative index of b in minorColumnKey:
704            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
705            poly signPoly = NULL;
706//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "poly signPoly = NULL;", k);
707            for (int r = 0; r < k; r++) {
708                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
709//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "MinorKey subMk = mk.getSubMinorKey(absoluteR, b);", k);
710                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
711                    PolyMinorValue mv;              // for storing all intermediate minors
712                    MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is mk with row absoluteR and column b omitted.
713//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "PolyMinorValue mv;", k);
714                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
715                        mv = cch.getValue(subMk);
716                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
717                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
718                                            // an impact on the internal ordering among cache entries.
719//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "cch.put(subMk, mv);", k);
720                    }
721                    else {
722                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
723                        // As this minor was not in the cache, we count the additions and
724                        // multiplications that we needed to do in the recursive call:
725                        m += mv.getMultiplications();
726                        s += mv.getAdditions();
727                    }
728                    // In any case, we count all nested operations in the accumulative counters:
729                    am += mv.getAccumulatedMultiplications();
730                    as += mv.getAccumulatedAdditions();
731                    pDelete(&signPoly);
732//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "pDelete(&signPoly);", k);
733                    signPoly = pISet(sign);
734//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "signPoly = pISet(sign);", k);
735                    poly temp = pp_Mult_qq(mv.getResult(), _polyMatrix[absoluteR][b], currRing);
736                    temp = p_Mult_q(signPoly, temp, currRing);
737                    result = p_Add_q(result, temp, currRing);
738//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d, case 2)", om_Info.UsedBytes, "ops(result, signPoly, mv.getResult(), _polyMatrix[absoluteR][b]);", k);
739                    signPoly = NULL;
740                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
741                }
742                sign = - sign; // alternating the sign
743            }
744        }
745        // Let's cache the newly computed minor:
746        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
747        s--; as--; // first addition was 0 + ..., so we do not count it
748        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
749        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
750        PolyMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
751//printf("\nMINOR to put: %s", pString(result));
752//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d)", om_Info.UsedBytes, "PolyMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);", k);
753        pDelete(&result); result = NULL;
754//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d)", om_Info.UsedBytes, "pDelete(&result);", k);
755        cch.put(mk, newMV); // Here's the actual put inside the cache.
756//omUpdateInfo(); printf("\nused bytes: %12ld, after: %s  (k = %d)", om_Info.UsedBytes, "cch.put(mk, newMV);", k);
757        return newMV;
758    }
759}
760
761#endif // HAVE_MINOR
762
763
764
765
Note: See TracBrowser for help on using the repository browser.