source: git/Singular/MinorProcessor.cc @ e2afced

spielwiese
Last change on this file since e2afced was e2afced, checked in by Frank Seelisch <seelisch@…>, 15 years ago
changes due to 32bit vs. 64bit (to make code run on compute servers) git-svn-id: file:///usr/local/Singular/svn/trunk@12170 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.2 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) {
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);
277}
278
279IntMinorValue IntMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices) {
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);
284}
285
286IntMinorValue IntMinorProcessor::getNextMinor() {
287    // computation without cache
288    return getMinorPrivate(_minorSize, _minor);
289}
290
291IntMinorValue IntMinorProcessor::getNextMinor(Cache<MinorKey, IntMinorValue>& c) {
292    // computation with cache
293    return getMinorPrivate(_minorSize, _minor, true, c);
294}
295
296IntMinorValue IntMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk) {
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);  // 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                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
327                }
328                sign = - sign; // alternating the sign
329            }
330        }
331        else {
332            b = - b - 1;
333            // This means that the best line is the column with absolute (0-based) index b.
334            // Using Laplace, the sign of the contributing minors must be iterating;
335            // the initial sign depends on the relative index of b in minorColumnKey:
336            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
337            for (int r = 0; r < k; r++) {
338                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
339                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
340                if (_matrix[absoluteR][b] != 0) { // Only then do we have to consider this sub-determinante.
341                    IntMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
342                    m += mv.getMultiplications();
343                    s += mv.getAdditions();
344                    am += mv.getAccumulatedMultiplications();
345                    as += mv.getAccumulatedAdditions();
346                    result += sign * mv.getResult() * _matrix[absoluteR][b]; // adding sub-determinante times matrix entry
347                                                                             // times appropriate sign
348                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
349                }
350                sign = - sign; // alternating the sign
351            }
352        }
353        s--; as--; // first addition was 0 + ..., so we do not count it
354        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
355        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
356        IntMinorValue newMV = IntMinorValue(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
357                                                                             // number of retrievals does not make sense, as we
358                                                                             // do not use a cache.
359        return newMV;
360    }
361}
362
363IntMinorValue IntMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk,
364                                                   const bool multipleMinors,
365                                                   Cache<MinorKey, IntMinorValue>& cch) {
366    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
367    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
368    if (k == 1) {
369        return IntMinorValue(_matrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
370                              0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
371    }
372    else {
373        int b = getBestLine(k, mk);                          // row or column with most zeros
374        int result = 0;                                      // This will contain the value of the minor.
375        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
376                                                             // ..."a*" for accumulated operation counters
377        IntMinorValue mv = IntMinorValue(0, 0, 0, 0, 0, 0, 0);     // for storing all intermediate minors
378        if (b >= 0) {
379            // This means that the best line is the row with absolute (0-based) index b.
380            // Using Laplace, the sign of the contributing minors must be iterating;
381            // the initial sign depends on the relative index of b in minorRowKey:
382            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
383            for (int c = 0; c < k; c++) {
384                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
385                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
386                if (_matrix[b][absoluteC] != 0) { // Only then do we have to consider this sub-determinante.
387                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
388                        mv = cch.getValue(subMk);
389                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
390                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
391                                            // an impact on the internal ordering among cache entries.
392                    }
393                    else {
394                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
395                        // As this minor was not in the cache, we count the additions and
396                        // multiplications that we needed to do in the recursive call:
397                        m += mv.getMultiplications();
398                        s += mv.getAdditions();
399                    }
400                    // In any case, we count all nested operations in the accumulative counters:
401                    am += mv.getAccumulatedMultiplications();
402                    as += mv.getAccumulatedAdditions();
403                    result += sign * mv.getResult() * _matrix[b][absoluteC]; // adding sub-determinante times matrix entry
404                                                                             // times appropriate sign
405                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
406                }
407                sign = - sign; // alternating the sign
408            }
409        }
410        else {
411            b = - b - 1;
412            // This means that the best line is the column with absolute (0-based) index b.
413            // Using Laplace, the sign of the contributing minors must be iterating;
414            // the initial sign depends on the relative index of b in minorColumnKey:
415            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
416            for (int r = 0; r < k; r++) {
417                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
418                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
419                if (_matrix[absoluteR][b] != 0) { // Only then do we have to consider this sub-determinante.
420                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
421                        mv = cch.getValue(subMk);
422                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
423                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
424                                            // an impact on the internal ordering among cache entries.
425                    }
426                    else {
427                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
428                        // As this minor was not in the cache, we count the additions and
429                        // multiplications that we needed to do in the recursive call:
430                        m += mv.getMultiplications();
431                        s += mv.getAdditions();
432                    }
433                    // In any case, we count all nested operations in the accumulative counters:
434                    am += mv.getAccumulatedMultiplications();
435                    as += mv.getAccumulatedAdditions();
436                    result += sign * mv.getResult() * _matrix[absoluteR][b]; // adding sub-determinante times matrix entry
437                                                                             // times appropriate sign
438                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
439                }
440                sign = - sign; // alternating the sign
441            }
442        }
443        // Let's cache the newly computed minor:
444        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
445        s--; as--; // first addition was 0 + ..., so we do not count it
446        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
447        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
448        IntMinorValue newMV = IntMinorValue(result, m, s, am, as, 1, potentialRetrievals);
449        cch.put(mk, newMV); // Here's the actual put inside the cache.
450        return newMV;
451    }
452}
453
454PolyMinorProcessor::PolyMinorProcessor () {
455    _polyMatrix = 0;
456}
457
458string PolyMinorProcessor::toString () const {
459    char h[32];
460    string t = "";
461    string s = "PolyMinorProcessor:";
462    s += "\n   matrix: ";
463    sprintf(h, "%d", _rows); s += h;
464    s += " x ";
465    sprintf(h, "%d", _columns); s += h;
466    int myIndexArray[500];
467    s += "\n   considered submatrix has row indices: ";
468    _container.getAbsoluteRowIndices(myIndexArray);
469    for (int k = 0; k < _containerRows; k++) {
470        if (k != 0) s += ", ";
471        sprintf(h, "%d", myIndexArray[k]); s += h;
472    }
473    s += " (first row of matrix has index 0)";
474    s += "\n   considered submatrix has column indices: ";
475    _container.getAbsoluteColumnIndices(myIndexArray);
476    for (int k = 0; k < _containerColumns; k++) {
477        if (k != 0) s += ", ";
478        sprintf(h, "%d", myIndexArray[k]); s += h;
479    }
480    s += " (first column of matrix has index 0)";
481    s += "\n   size of considered minor(s): ";
482    sprintf(h, "%d", _minorSize); s += h;
483    s += "x";
484    s += h;
485    return s;
486}
487
488bool PolyMinorProcessor::isEntryZero (const int absoluteRowIndex, const int absoluteColumnIndex) const
489{
490  poly zeroPoly = pISet(0);
491  return pEqualPolys(_polyMatrix[absoluteRowIndex][absoluteColumnIndex], zeroPoly);
492}
493
494PolyMinorProcessor::~PolyMinorProcessor() {
495    // free memory of _polyMatrix
496    for (int i = 0; i < _rows; i++) {
497        for (int j = 0; j < _columns; j++) {
498            p_Delete(&_polyMatrix[i][j], currRing);
499        }
500        delete [] _polyMatrix[i]; _polyMatrix[i] = 0;
501    }
502    delete [] _polyMatrix; _polyMatrix = 0;
503}
504
505void PolyMinorProcessor::defineMatrix (const int numberOfRows, const int numberOfColumns, const poly* polyMatrix) {
506    // free memory of _polyMatrix
507    for (int i = 0; i < _rows; i++) {
508        for (int j = 0; j < _columns; j++)
509            p_Delete(&_polyMatrix[i][j], currRing);
510        delete [] _polyMatrix[i]; _polyMatrix[i] = 0;
511    }
512    delete [] _polyMatrix; _polyMatrix = 0;
513
514    _rows = numberOfRows;
515    _columns = numberOfColumns;
516
517    // allocate memory for new entries in _matrix
518    _polyMatrix = new poly*[_rows];
519    for (int i = 0; i < _rows; i++) _polyMatrix[i] = new poly[_columns];
520
521    // copying values from one-dimensional method parameter "matrix"
522    for (int r = 0; r < _rows; r++)
523        for (int c = 0; c < _columns; c++)
524            _polyMatrix[r][c] = pCopy(polyMatrix[r * _columns + c]);
525}
526
527PolyMinorValue PolyMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices,
528                                            Cache<MinorKey, PolyMinorValue>& c) {
529    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
530    _minorSize = dimension;
531    // call a helper method which recursively computes the minor using the cache c:
532    return getMinorPrivate(dimension, _container, false, c);
533}
534
535PolyMinorValue PolyMinorProcessor::getMinor(const int dimension, const int* rowIndices, const int* columnIndices) {
536    defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
537    _minorSize = dimension;
538    // call a helper method which recursively computes the minor (without using a cache):
539    return getMinorPrivate(_minorSize, _container);
540}
541
542PolyMinorValue PolyMinorProcessor::getNextMinor() {
543    // computation without cache
544    return getMinorPrivate(_minorSize, _minor);
545}
546
547PolyMinorValue PolyMinorProcessor::getNextMinor(Cache<MinorKey, PolyMinorValue>& c) {
548    // computation with cache
549    return getMinorPrivate(_minorSize, _minor, true, c);
550}
551
552// performs the assignment a = a + (b * c * d)
553// c and d must neither be destroyed nor modified
554// b may be destroyed
555// a must finally contain the new value
556poly ops(poly a, poly b, poly c, poly d)
557{
558  return p_Add_q(a, p_Mult_q(b, pp_Mult_qq(c, d, currRing), currRing), currRing);
559}
560
561PolyMinorValue PolyMinorProcessor::getMinorPrivate(const int k, const MinorKey& mk) {
562    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
563    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
564    if (k == 1) {
565        return PolyMinorValue(_polyMatrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
566                              0, 0, 0, 0, -1, -1); // "-1" is to signal that any statistics about the
567                                                   // number of retrievals does not make sense, as we do not use a cache.
568    }
569    else {
570        // Here, the minor must be 2x2 or larger.
571        int b = getBestLine(k, mk);                          // row or column with most zeros
572        poly result = pISet(0);                              // This will contain the value of the minor.
573        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
574                                                             // ..."a*" for accumulated operation counters
575        if (b >= 0) {
576            // This means that the best line is the row with absolute (0-based) index b.
577            // Using Laplace, the sign of the contributing minors must be iterating;
578            // the initial sign depends on the relative index of b in minorRowKey:
579            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
580            poly signPoly = pISet(sign);
581            for (int c = 0; c < k; c++) {
582                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
583                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
584                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
585                    PolyMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
586                    m += mv.getMultiplications();
587                    s += mv.getAdditions();
588                    am += mv.getAccumulatedMultiplications();
589                    as += mv.getAccumulatedAdditions();
590                    result = ops(result, signPoly, mv.getResult(), _polyMatrix[b][absoluteC]); // adding sub-determinante times matrix entry
591                                                                                               // times appropriate sign
592                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
593                }
594                sign = - sign; // alternating the sign
595                signPoly = pISet(sign);
596            }
597            // p_Delete(&signPoly, currRing);
598        }
599        else {
600            b = - b - 1;
601            // This means that the best line is the column with absolute (0-based) index b.
602            // Using Laplace, the sign of the contributing minors must be iterating;
603            // the initial sign depends on the relative index of b in minorColumnKey:
604            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
605            poly signPoly = pISet(sign);
606            for (int r = 0; r < k; r++) {
607                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
608                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
609                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
610                    PolyMinorValue mv = getMinorPrivate(k - 1, subMk);  // recursive call
611                    m += mv.getMultiplications();
612                    s += mv.getAdditions();
613                    am += mv.getAccumulatedMultiplications();
614                    as += mv.getAccumulatedAdditions();
615                    result = ops(result, signPoly, mv.getResult(), _polyMatrix[absoluteR][b]); // adding sub-determinante times matrix entry
616                                                                                               // times appropriate sign
617                    s++; m++; as++, am++; // This is for the addition and multiplication in the previous line of code.
618                }
619                sign = - sign; // alternating the sign
620                signPoly = pISet(sign);
621            }
622            // p_Delete(&signPoly, currRing);
623        }
624        s--; as--; // first addition was 0 + ..., so we do not count it
625        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
626        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
627        PolyMinorValue newMV(result, m, s, am, as, -1, -1); // "-1" is to signal that any statistics about the
628                                                            // number of retrievals does not make sense, as we
629                                                            // do not use a cache.
630        // p_Delete(&result, currRing);
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    assert(k > 0); // k is the minor's dimension; the minor must be at least 1x1
639    // The method works by recursion, and using Lapace's Theorem along the row/column with the most zeros.
640    if (k == 1) {
641        return PolyMinorValue(_polyMatrix[mk.getAbsoluteRowIndex(0)][mk.getAbsoluteColumnIndex(0)],
642                              0, 0, 0, 0, -1, -1); // we set "-1" as, for k == 1, we do not have any cache retrievals
643    }
644    else {
645        int b = getBestLine(k, mk);                          // row or column with most zeros
646        poly result = pISet(0);                              // This will contain the value of the minor.
647        int s = 0; int m = 0; int as = 0; int am = 0;        // counters for additions and multiplications,
648                                                             // ..."a*" for accumulated operation counters
649        PolyMinorValue mv = PolyMinorValue(0, 0, 0, 0, 0, 0, 0);     // for storing all intermediate minors
650        if (b >= 0) {
651            // This means that the best line is the row with absolute (0-based) index b.
652            // Using Laplace, the sign of the contributing minors must be iterating;
653            // the initial sign depends on the relative index of b in minorRowKey:
654            int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
655            poly signPoly = pISet(sign);
656            for (int c = 0; c < k; c++) {
657                int absoluteC = mk.getAbsoluteColumnIndex(c);      // This iterates over all involved columns.
658                MinorKey subMk = mk.getSubMinorKey(b, absoluteC);  // This is MinorKey when we omit row b and column absoluteC.
659                if (!isEntryZero(b, absoluteC)) { // Only then do we have to consider this sub-determinante.
660                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
661                        mv = cch.getValue(subMk);
662                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
663                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
664                                            // an impact on the internal ordering among cache entries.
665                    }
666                    else {
667                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
668                        // As this minor was not in the cache, we count the additions and
669                        // multiplications that we needed to do in the recursive call:
670                        m += mv.getMultiplications();
671                        s += mv.getAdditions();
672                    }
673                    // In any case, we count all nested operations in the accumulative counters:
674                    am += mv.getAccumulatedMultiplications();
675                    as += mv.getAccumulatedAdditions();
676                    result = ops(result, signPoly, mv.getResult(), _polyMatrix[b][absoluteC]); // adding sub-determinante times matrix entry
677                                                                                               // times appropriate sign
678                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
679                }
680                sign = - sign; // alternating the sign
681                signPoly = pISet(sign);
682            }
683        }
684        else {
685            b = - b - 1;
686            // This means that the best line is the column with absolute (0-based) index b.
687            // Using Laplace, the sign of the contributing minors must be iterating;
688            // the initial sign depends on the relative index of b in minorColumnKey:
689            int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
690            poly signPoly = pISet(sign);
691            for (int r = 0; r < k; r++) {
692                int absoluteR = mk.getAbsoluteRowIndex(r);        // This iterates over all involved rows.
693                MinorKey subMk = mk.getSubMinorKey(absoluteR, b); // This is MinorKey when we omit row absoluteR and column b.
694                if (!isEntryZero(absoluteR, b)) { // Only then do we have to consider this sub-determinante.
695                    if (cch.hasKey(subMk)) { // trying to find the result in the cache
696                        mv = cch.getValue(subMk);
697                        mv.incrementRetrievals(); // once more, we made use of the cached value for key mk
698                        cch.put(subMk, mv); // We need to do this with "put", as the (altered) number of retrievals may have
699                                            // an impact on the internal ordering among cache entries.
700                    }
701                    else {
702                        mv = getMinorPrivate(k - 1, subMk, multipleMinors, cch); // recursive call
703                        // As this minor was not in the cache, we count the additions and
704                        // multiplications that we needed to do in the recursive call:
705                        m += mv.getMultiplications();
706                        s += mv.getAdditions();
707                    }
708                    // In any case, we count all nested operations in the accumulative counters:
709                    am += mv.getAccumulatedMultiplications();
710                    as += mv.getAccumulatedAdditions();
711                    result = ops(result, signPoly, mv.getResult(), _polyMatrix[absoluteR][b]); // adding sub-determinante times matrix entry
712                                                                                               // times appropriate sign
713                    s++; m++; as++; am++; // This is for the addition and multiplication in the previous line of code.
714                }
715                sign = - sign; // alternating the sign
716                signPoly = pISet(sign);
717            }
718        }
719        // Let's cache the newly computed minor:
720        int potentialRetrievals = NumberOfRetrievals(_containerRows, _containerColumns, _minorSize, k, multipleMinors);
721        s--; as--; // first addition was 0 + ..., so we do not count it
722        if (s < 0) s = 0; // may happen when all subminors are zero and no addition needs to be performed
723        if (as < 0) as = 0; // may happen when all subminors are zero and no addition needs to be performed
724        PolyMinorValue newMV = PolyMinorValue(result, m, s, am, as, 1, potentialRetrievals);
725        cch.put(mk, newMV); // Here's the actual put inside the cache.
726        return newMV;
727    }
728}
729
730#endif // HAVE_MINOR
731
732
733
734
Note: See TracBrowser for help on using the repository browser.