# source:git/Singular/MinorProcessor.cc@40349c

spielwiese
Last change on this file since 40349c was 40349c, checked in by Frank Seelisch <seelisch@…>, 14 years ago
dealing with polynomials / pDelete git-svn-id: file:///usr/local/Singular/svn/trunk@12210 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 37.9 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                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
318                if (_matrix[b][absoluteC] != 0) { // Only then do we have to consider this sub-determinante.
319                    IntMinorValue mv = getMinorPrivate(k - 1, subMk, characteristic);  // recursive call
320                    m += mv.getMultiplications();
322                    am += mv.getAccumulatedMultiplications();
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                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
341                if (_matrix[absoluteR][b] != 0) { // Only then do we have to consider this sub-determinante.
342                    IntMinorValue mv = getMinorPrivate(k - 1, subMk, characteristic);  // recursive call
343                    m += mv.getMultiplications();
345                    am += mv.getAccumulatedMultiplications();
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                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
389                if (_matrix[b][absoluteC] != 0) { // Only then do we have to consider this sub-determinante.
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();
402                    }
403                    // In any case, we count all nested operations in the accumulative counters:
404                    am += mv.getAccumulatedMultiplications();
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                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
423                if (_matrix[absoluteR][b] != 0) { // Only then do we have to consider this sub-determinante.
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();
436                    }
437                    // In any case, we count all nested operations in the accumulative counters:
438                    am += mv.getAccumulatedMultiplications();
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
556// performs the assignment a = a + (b * c * d)
557// c and d must neither be destroyed nor modified
558// the old value of a and b will be destroyed
559// a will finally contain the new value
560void ops(poly a, poly b, poly c, poly d)
561{
562  a = p_Add_q(a, p_Mult_q(b, pp_Mult_qq(c, d, currRing), currRing), currRing);
563}
564
565PolyMinorValue PolyMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk) {
566    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
567    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
568    if (k == 1) {
569        return PolyMinorValue(_polyMatrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
570                              0, 0, 0, 0, -1, -1); // "-1" is to signal that any statistics about the
571                                                   // number of retrievals does not make sense, as we do not use a cache.
572    }
573    else {
574        // Here, the minor must be 2x2 or larger.
575        int b = getBestLine(k, mk);                          // row or column with most zeros
576        poly result = pISet(0);                              // This will contain the value of the minor.
577        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
578                                                             // ..."a*" for accumulated operation counters
579        if (b >= 0) {
580            // This means that the best line is the row with absolute (0-based) index b.
581            // Using Laplace, the sign of the contributing minors must be iterating;
582            // the initial sign depends on the relative index of b in minorRowKey:
583            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
584            poly signPoly = NULL;
585            for (int c = 0; c < k; c++) {
586                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
587                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
588                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
589                    PolyMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
590                    m += mv.getMultiplications();
592                    am += mv.getAccumulatedMultiplications();
594                    pDelete(&signPoly);
595                    signPoly = pISet(sign);
596                    ops(result, signPoly, mv.getResult(), _polyMatrix[b][absoluteC]); // adding in "result" the product of sign,
597                                                                                      // sub-determinante, and matrix entry
598                    signPoly = NULL;
599                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
600                }
601                sign = - sign; // alternating the sign
602            }
603        }
604        else {
605            b = - b - 1;
606            // This means that the best line is the column with absolute (0-based) index b.
607            // Using Laplace, the sign of the contributing minors must be iterating;
608            // the initial sign depends on the relative index of b in minorColumnKey:
609            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
610            poly signPoly = NULL;
611            for (int r = 0; r < k; r++) {
612                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
613                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
614                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
615                    PolyMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
616                    m += mv.getMultiplications();
618                    am += mv.getAccumulatedMultiplications();
620                    pDelete(&signPoly);
621                    signPoly = pISet(sign);
622                    ops(result, signPoly, mv.getResult(), _polyMatrix[absoluteR][b]); // adding in "result" the product of sign,
623                                                                                      // sub-determinante, and matrix entry
624                    signPoly = NULL;
625                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
626                }
627                sign = - sign; // alternating the sign
628            }
629        }
630        s--; as--; // first addition was 0 + ..., so we do not count it
631        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
632        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
633        PolyMinorValue newMV(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
634                                                            // number of retrievals does not make sense, as we
635                                                            // do not use a cache.
636        pDelete(&result);
637        return newMV;
638    }
639}
640
641PolyMinorValue PolyMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk,
642                                                   const bool multipleMinors,
643                                                   Cache<MinorKey, PolyMinorValue>& cch) {
644    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
645    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
646    if (k == 1) {
647        return PolyMinorValue(_polyMatrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
648                              0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
649    }
650    else {
651        int b = getBestLine(k, mk);                          // row or column with most zeros
652        poly result = pISet(0);                              // This will contain the value of the minor.
653        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
654                                                             // ..."a*" for accumulated operation counters
655        if (b >= 0) {
656            // This means that the best line is the row with absolute (0-based) index b.
657            // Using Laplace, the sign of the contributing minors must be iterating;
658            // the initial sign depends on the relative index of b in minorRowKey:
659            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
660            poly signPoly = NULL;
661            for (int c = 0; c < k; c++) {
662                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
663                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
664                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
665                    PolyMinorValue mv;              // for storing all intermediate minors
666                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
667                        mv = cch.getValue(subMk);
668                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
669                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
670                                            // an impact on the internal ordering among cache entries.
671                    }
672                    else {
673                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
674                        // As this minor was not in the cache, we count the additions and
675                        // multiplications that we needed to do in the recursive call:
676                        m += mv.getMultiplications();
678                    }
679                    // In any case, we count all nested operations in the accumulative counters:
680                    am += mv.getAccumulatedMultiplications();
682                    pDelete(&signPoly);
683                    signPoly = pISet(sign);
684                    ops(result, signPoly, mv.getResult(), _polyMatrix[b][absoluteC]); // adding in "result" the product of sign,
685                                                                                      // sub-determinante, and matrix entry
686                    signPoly = NULL;
687                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
688                }
689                sign = - sign; // alternating the sign
690            }
691        }
692        else {
693            b = - b - 1;
694            // This means that the best line is the column with absolute (0-based) index b.
695            // Using Laplace, the sign of the contributing minors must be iterating;
696            // the initial sign depends on the relative index of b in minorColumnKey:
697            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
698            poly signPoly = NULL;
699            for (int r = 0; r < k; r++) {
700                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
701                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
702                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
703                    PolyMinorValue mv;              // for storing all intermediate minors
704                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
705                        mv = cch.getValue(subMk);
706                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
707                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
708                                            // an impact on the internal ordering among cache entries.
709                    }
710                    else {
711                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
712                        // As this minor was not in the cache, we count the additions and
713                        // multiplications that we needed to do in the recursive call:
714                        m += mv.getMultiplications();
716                    }
717                    // In any case, we count all nested operations in the accumulative counters:
718                    am += mv.getAccumulatedMultiplications();
720                    pDelete(&signPoly);
721                    signPoly = pISet(sign);
722                    ops(result, signPoly, mv.getResult(), _polyMatrix[absoluteR][b]); // adding in "result" the product of sign,
723                                                                                      // sub-determinante, and matrix entry
724                    signPoly = NULL;
725                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
726                }
727                sign = - sign; // alternating the sign
728            }
729        }
730        // Let's cache the newly computed minor:
731        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
732        s--; as--; // first addition was 0 + ..., so we do not count it
733        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
734        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
735        PolyMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
736        pDelete(&result);
737        cch.put(mk, newMV); // Here's the actual put inside the cache.
738        return newMV;
739    }
740}
741
742#endif // HAVE_MINOR
743
744
745
746
Note: See TracBrowser for help on using the repository browser.