My Project
Loading...
Searching...
No Matches
MinorProcessor.cc
Go to the documentation of this file.
1
2
3
4#include "kernel/mod2.h"
5
7
8#include "polys/kbuckets.h"
9
10#include "kernel/structs.h"
11#include "kernel/polys.h"
13
14#include "kernel/ideals.h"
15
16using namespace std;
17
18#ifdef COUNT_AND_PRINT_OPERATIONS
19VAR long addsPoly = 0; /* for the number of additions of two polynomials */
20VAR long multsPoly = 0; /* for the number of multiplications of two polynomials */
21VAR long addsPolyForDiv = 0; /* for the number of additions of two polynomials for
22 polynomial division part */
23VAR long multsPolyForDiv = 0; /* for the number of multiplications of two polynomials
24 for polynomial division part */
25VAR long multsMon = 0; /* for the number of multiplications of two monomials */
26VAR long multsMonForDiv = 0; /* for the number of m-m-multiplications for polynomial
27 division part */
28VAR long savedMultsMFD = 0; /* number of m-m-multiplications that could be saved
29 when polynomial division would be optimal
30 (if p / t1 = t2 + ..., then t1 * t2 = LT(p), i.e.,
31 this multiplication need not be performed which
32 would save one m-m-multiplication) */
33VAR long divsMon = 0; /* for the number of divisions of two monomials;
34 these are all guaranteed to work, i.e., m1/m2 only
35 when exponentVector(m1) >= exponentVector(m2) */
36void printCounters (char* prefix, bool resetToZero)
37{
38 printf("%s [p+p(div) | p*p(div) | m*m(div, -save) | m/m ]", prefix);
39 printf(" = [%ld(%ld) | %ld(%ld) | %ld(%d, -%ld) | %ld]\n",
40 addsPoly, addsPolyForDiv, multsPoly, multsPolyForDiv,
41 multsMon, multsMonForDiv, savedMultsMFD, divsMon);
42 if (resetToZero)
43 {
44 multsMon = 0; addsPoly = 0; multsPoly = 0; divsMon = 0;
45 savedMultsMFD = 0; multsMonForDiv = 0; addsPolyForDiv = 0;
46 multsPolyForDiv = 0;
47 }
48}
49#endif
50/* COUNT_AND_PRINT_OPERATIONS */
51
53{
54 PrintS(this->toString().c_str());
55}
56
57int MinorProcessor::getBestLine (const int k, const MinorKey& mk) const
58{
59 /* This method identifies the row or column with the most zeros.
60 The returned index (bestIndex) is absolute within the pre-
61 defined matrix.
62 If some row has the most zeros, then the absolute (0-based)
63 row index is returned.
64 If, contrariwise, some column has the most zeros, then -1 minus
65 the absolute (0-based) column index is returned. */
66 int numberOfZeros = 0;
67 int bestIndex = 100000; /* We start with an invalid row/column index. */
68 int maxNumberOfZeros = -1; /* We update this variable whenever we find
69 a new so-far optimal row or column. */
70 for (int r = 0; r < k; r++)
71 {
72 /* iterate through all k rows of the momentary minor */
73 int absoluteR = mk.getAbsoluteRowIndex(r);
74 numberOfZeros = 0;
75 for (int c = 0; c < k; c++)
76 {
77 int absoluteC = mk.getAbsoluteColumnIndex(c);
78 if (isEntryZero(absoluteR, absoluteC)) numberOfZeros++;
79 }
80 if (numberOfZeros > maxNumberOfZeros)
81 {
82 /* We found a new best line which is a row. */
83 bestIndex = absoluteR;
84 maxNumberOfZeros = numberOfZeros;
85 }
86 };
87 for (int c = 0; c < k; c++)
88 {
89 int absoluteC = mk.getAbsoluteColumnIndex(c);
90 numberOfZeros = 0;
91 for (int r = 0; r < k; r++)
92 {
93 int absoluteR = mk.getAbsoluteRowIndex(r);
94 if (isEntryZero(absoluteR, absoluteC)) numberOfZeros++;
95 }
96 if (numberOfZeros > maxNumberOfZeros)
97 {
98 /* We found a new best line which is a column. So we transform
99 the return value. Note that we can easily retrieve absoluteC
100 from bestLine: absoluteC = - 1 - bestLine. */
101 bestIndex = - absoluteC - 1;
102 maxNumberOfZeros = numberOfZeros;
103 }
104 };
105 return bestIndex;
106}
107
108void MinorProcessor::setMinorSize(const int minorSize)
109{
110 _minorSize = minorSize;
111 _minor.reset();
112}
113
115{
116 return setNextKeys(_minorSize);
117}
118
119void MinorProcessor::getCurrentRowIndices(int* const target) const
120{
121 return _minor.getAbsoluteRowIndices(target);
122}
123
124void MinorProcessor::getCurrentColumnIndices(int* const target) const
125{
126 return _minor.getAbsoluteColumnIndices(target);
127}
128
129void MinorProcessor::defineSubMatrix(const int numberOfRows,
130 const int* rowIndices,
131 const int numberOfColumns,
132 const int* columnIndices)
133{
134 /* The method assumes ascending row and column indices in the
135 two argument arrays. These indices are understood to be zero-based.
136 The method will set the two arrays of ints in _container.
137 Example: The indices 0, 2, 3, 7 will be converted to an array with
138 one int representing the binary number 10001101
139 (check bits from right to left). */
140
141 _containerRows = numberOfRows;
142 int highestRowIndex = rowIndices[numberOfRows - 1];
143 int rowBlockCount = (highestRowIndex / 32) + 1;
144 unsigned *rowBlocks=(unsigned*)omAlloc(rowBlockCount*sizeof(unsigned));
145 for (int i = 0; i < rowBlockCount; i++) rowBlocks[i] = 0;
146 for (int i = 0; i < numberOfRows; i++)
147 {
148 int blockIndex = rowIndices[i] / 32;
149 int offset = rowIndices[i] % 32;
150 rowBlocks[blockIndex] += (1 << offset);
151 }
152
153 _containerColumns = numberOfColumns;
154 int highestColumnIndex = columnIndices[numberOfColumns - 1];
155 int columnBlockCount = (highestColumnIndex / 32) + 1;
156 unsigned *columnBlocks=(unsigned*)omAlloc0(columnBlockCount*sizeof(unsigned));
157 for (int i = 0; i < numberOfColumns; i++)
158 {
159 int blockIndex = columnIndices[i] / 32;
160 int offset = columnIndices[i] % 32;
161 columnBlocks[blockIndex] += (1 << offset);
162 }
163
164 _container.set(rowBlockCount, rowBlocks, columnBlockCount, columnBlocks);
165 omFree(columnBlocks);
166 omFree(rowBlocks);
167}
168
170{
171 /* This method moves _minor to the next valid (k x k)-minor within
172 _container. It returns true iff this is successful, i.e. iff
173 _minor did not already encode the terminal (k x k)-minor. */
174 if (_minor.compare(MinorKey(0, 0, 0, 0)) == 0)
175 {
176 /* This means that we haven't started yet. Thus, we are about
177 to compute the first (k x k)-minor. */
180 return true;
181 }
183 {
184 /* Here we were able to pick a next subset of columns
185 within the same subset of rows. */
186 return true;
187 }
189 {
190 /* Here we were not able to pick a next subset of columns
191 within the same subset of rows. But we could pick a next
192 subset of rows. We must hence reset the subset of columns: */
194 return true;
195 }
196 else
197 {
198 /* We were neither able to pick a next subset
199 of columns nor of rows. I.e., we have iterated through
200 all sensible choices of subsets of rows and columns. */
201 return false;
202 }
203}
204
205bool MinorProcessor::isEntryZero (const int /*absoluteRowIndex*/,
206 const int /*absoluteColumnIndex*/) const
207{
208 assume(false);
209 return false;
210}
211
213{
214 assume(false);
215 return "";
216}
217
218int MinorProcessor::IOverJ(const int i, const int j)
219{
220 /* This is a non-recursive implementation. */
221 assume( (i >= 0) && (j >= 0) && (i >= j));
222 if (j == 0 || i == j) return 1;
223 #if 0
224 int result = 1;
225 for (int k = i - j + 1; k <= i; k++) result *= k;
226 /* Now, result = (i - j + 1) * ... * i. */
227 for (int k = 2; k <= j; k++) result /= k;
228 /* Now, result = (i - j + 1) * ... * i / 1 / 2 ...
229 ... / j = i! / j! / (i - j)!. */
230 return result;
231 #endif
232 return binom(i,j);
233}
234
236{
237 /* This is a non-recursive implementation. */
238 assume(i >= 0);
239 int result = 1;
240 for (int j = 1; j <= i; j++) result *= j;
241 // Now, result = 1 * 2 * ... * i = i!
242 return result;
243}
244
245int MinorProcessor::NumberOfRetrievals (const int rows, const int columns,
246 const int containerMinorSize,
247 const int minorSize,
248 const bool multipleMinors)
249{
250 /* This method computes the number of potential retrievals
251 of a single minor when computing all minors of a given size
252 within a matrix of given size. */
253 int result = 0;
254 if (multipleMinors)
255 {
256 /* Here, we would like to compute all minors of size
257 containerMinorSize x containerMinorSize in a matrix
258 of size rows x columns.
259 Then, we need to retrieve any minor of size
260 minorSize x minorSize exactly n times, where n is as
261 follows: */
262 result = IOverJ(rows - minorSize, containerMinorSize - minorSize)
263 * IOverJ(columns - minorSize, containerMinorSize - minorSize)
264 * Faculty(containerMinorSize - minorSize);
265 }
266 else
267 {
268 /* Here, we would like to compute just one minor of size
269 containerMinorSize x containerMinorSize. Then, we need
270 to retrieve any minor of size minorSize x minorSize exactly
271 (containerMinorSize - minorSize)! times: */
272 result = Faculty(containerMinorSize - minorSize);
273 }
274 return result;
275}
276
278 _container(0, NULL, 0, NULL),
279 _minor(0, NULL, 0, NULL),
280 _containerRows(0),
281 _containerColumns(0),
282 _minorSize(0),
283 _rows(0),
284 _columns(0)
285{
286}
287
289
291{
292 _intMatrix = 0;
293}
294
296{
297 char h[32];
298 string t = "";
299 string s = "IntMinorProcessor:";
300 s += "\n matrix: ";
301 sprintf(h, "%d", _rows); s += h;
302 s += " x ";
303 sprintf(h, "%d", _columns); s += h;
304 for (int r = 0; r < _rows; r++)
305 {
306 s += "\n ";
307 for (int c = 0; c < _columns; c++)
308 {
309 sprintf(h, "%d", getEntry(r, c)); t = h;
310 for (int k = 0; k < int(4 - strlen(h)); k++) s += " ";
311 s += t;
312 }
313 }
314 int myIndexArray[500];
315 s += "\n considered submatrix has row indices: ";
316 _container.getAbsoluteRowIndices(myIndexArray);
317 for (int k = 0; k < _containerRows; k++)
318 {
319 if (k != 0) s += ", ";
320 sprintf(h, "%d", myIndexArray[k]); s += h;
321 }
322 s += " (first row of matrix has index 0)";
323 s += "\n considered submatrix has column indices: ";
325 for (int k = 0; k < _containerColumns; k++)
326 {
327 if (k != 0) s += ", ";
328 sprintf(h, "%d", myIndexArray[k]); s += h;
329 }
330 s += " (first column of matrix has index 0)";
331 s += "\n size of considered minor(s): ";
332 sprintf(h, "%d", _minorSize); s += h;
333 s += "x";
334 s += h;
335 return s;
336}
337
339{
340 /* free memory of _intMatrix */
341 delete [] _intMatrix; _intMatrix = 0;
342}
343
344bool IntMinorProcessor::isEntryZero (const int absoluteRowIndex,
345 const int absoluteColumnIndex) const
346{
347 return getEntry(absoluteRowIndex, absoluteColumnIndex) == 0;
348}
349
350void IntMinorProcessor::defineMatrix (const int numberOfRows,
351 const int numberOfColumns,
352 const int* matrix)
353{
354 /* free memory of _intMatrix */
356
357 _rows = numberOfRows;
358 _columns = numberOfColumns;
359
360 /* allocate memory for new entries in _intMatrix */
361 int n = _rows * _columns;
362 _intMatrix = (int*)omAlloc(n*sizeof(int));
363
364 /* copying values from one-dimensional method
365 parameter "matrix" */
366 for (int i = 0; i < n; i++)
367 _intMatrix[i] = matrix[i];
368}
369
371 const int* rowIndices,
372 const int* columnIndices,
374 const int characteristic,
375 const ideal& iSB)
376{
377 defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
379 /* call a helper method which recursively computes the minor using the
380 cache c: */
382 characteristic, iSB);
383}
384
386 const int* rowIndices,
387 const int* columnIndices,
388 const int characteristic,
389 const ideal& iSB,
390 const char* algorithm)
391{
392 defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
394
395 /* call a helper method which computes the minor (without a cache): */
396 if (strcmp(algorithm, "Laplace") == 0)
397 return getMinorPrivateLaplace(_minorSize, _container, characteristic,
398 iSB);
399 else if (strcmp(algorithm, "Bareiss") == 0)
400 return getMinorPrivateBareiss(_minorSize, _container, characteristic,
401 iSB);
402 else assume(false);
403
404 /* The following code is never reached and just there to make the
405 compiler happy: */
406 return IntMinorValue();
407}
408
410 const ideal& iSB,
411 const char* algorithm)
412{
413 /* call a helper method which computes the minor (without a cache): */
414 if (strcmp(algorithm, "Laplace") == 0)
415 return getMinorPrivateLaplace(_minorSize, _minor, characteristic, iSB);
416 else if (strcmp(algorithm, "Bareiss") == 0)
417 return getMinorPrivateBareiss(_minorSize, _minor, characteristic, iSB);
418 else assume(false);
419
420 /* The following code is never reached and just there to make the
421 compiler happy: */
422 return IntMinorValue();
423}
424
426 IntMinorValue>& c,
427 const int characteristic,
428 const ideal& iSB)
429{
430 /* computation with cache */
431 return getMinorPrivateLaplace(_minorSize, _minor, true, c, characteristic,
432 iSB);
433}
434
435/* computes the reduction of an integer i modulo an ideal
436 which captures a std basis */
437int getReduction (const int i, const ideal& iSB)
438{
439 if (i == 0) return 0;
440 poly f = pISet(i);
441 poly g = kNF(iSB, currRing->qideal, f);
442 int result = 0;
443 if (g != NULL) result = n_Int(pGetCoeff(g), currRing->cf);
444 pDelete(&f);
445 pDelete(&g);
446 return result;
447}
448
450 const int k,
451 const MinorKey& mk,
452 const int characteristic,
453 const ideal& iSB)
454{
455 assume(k > 0); /* k is the minor's dimension; the minor must be at least
456 1x1 */
457 /* The method works by recursion, and using Lapace's Theorem along the
458 row/column with the most zeros. */
459 if (k == 1)
460 {
462 if (characteristic != 0) e = e % characteristic;
463 if (iSB != 0) e = getReduction(e, iSB);
464 return IntMinorValue(e, 0, 0, 0, 0, -1, -1); /* "-1" is to signal that any
465 statistics about the number
466 of retrievals does not make
467 sense, as we do not use a
468 cache. */
469 }
470 else
471 {
472 /* Here, the minor must be 2x2 or larger. */
473 int b = getBestLine(k, mk); /* row or column with most
474 zeros */
475 int result = 0; /* This will contain the
476 value of the minor. */
477 int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions and
478 multiplications, ..."a*"
479 for accumulated operation
480 counters */
481 bool hadNonZeroEntry = false;
482 if (b >= 0)
483 {
484 /* This means that the best line is the row with absolute (0-based)
485 index b.
486 Using Laplace, the sign of the contributing minors must be iterating;
487 the initial sign depends on the relative index of b in minorRowKey: */
488 int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
489 for (int c = 0; c < k; c++) /* This iterates over all involved columns. */
490 {
491 int absoluteC = mk.getAbsoluteColumnIndex(c);
492 if (getEntry(b, absoluteC) != 0) /* Only then do we have to consider
493 this sub-determinante. */
494 {
495 hadNonZeroEntry = true;
496 /* Next MinorKey is mk with row b and column absoluteC omitted: */
497 MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
498 IntMinorValue mv = getMinorPrivateLaplace(k - 1, subMk,
499 characteristic, iSB); /* recursive call */
500 m += mv.getMultiplications();
501 s += mv.getAdditions();
503 as += mv.getAccumulatedAdditions();
504 /* adding sub-determinante times matrix entry
505 times appropriate sign: */
506 result += sign * mv.getResult() * getEntry(b, absoluteC);
507
508 if (characteristic != 0) result = result % characteristic;
509 s++; m++; as++, am++; /* This is for the last addition and
510 multiplication. */
511 }
512 sign = - sign; /* alternating the sign */
513 }
514 }
515 else
516 {
517 b = - b - 1;
518 /* This means that the best line is the column with absolute (0-based)
519 index b.
520 Using Laplace, the sign of the contributing minors must be iterating;
521 the initial sign depends on the relative index of b in
522 minorColumnKey: */
523 int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
524 for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
525 {
526 int absoluteR = mk.getAbsoluteRowIndex(r);
527 if (getEntry(absoluteR, b) != 0) /* Only then do we have to consider
528 this sub-determinante. */
529 {
530 hadNonZeroEntry = true;
531 /* Next MinorKey is mk with row absoluteR and column b omitted. */
532 MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
533 IntMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, characteristic, iSB); /* recursive call */
534 m += mv.getMultiplications();
535 s += mv.getAdditions();
537 as += mv.getAccumulatedAdditions();
538 /* adding sub-determinante times matrix entry
539 times appropriate sign: */
540 result += sign * mv.getResult() * getEntry(absoluteR, b);
541 if (characteristic != 0) result = result % characteristic;
542 s++; m++; as++, am++; /* This is for the last addition and
543 multiplication. */
544 }
545 sign = - sign; /* alternating the sign */
546 }
547 }
548 if (hadNonZeroEntry)
549 {
550 s--; as--; /* first addition was 0 + ..., so we do not count it */
551 }
552 if (s < 0) s = 0; /* may happen when all subminors are zero and no
553 addition needs to be performed */
554 if (as < 0) as = 0; /* may happen when all subminors are zero and no
555 addition needs to be performed */
556 if (iSB != 0) result = getReduction(result, iSB);
557 IntMinorValue newMV(result, m, s, am, as, -1, -1);
558 /* "-1" is to signal that any statistics about the number of retrievals
559 does not make sense, as we do not use a cache. */
560 return newMV;
561 }
562}
563
564/* This method can only be used in the case of coefficients
565 coming from a field or at least from an integral domain. */
567 const int k,
568 const MinorKey& mk,
569 const int characteristic,
570 const ideal& iSB)
571{
572 assume(k > 0); /* k is the minor's dimension; the minor must be at least
573 1x1 */
574 int *theRows=(int*)omAlloc(k*sizeof(int));
575 mk.getAbsoluteRowIndices(theRows);
576 int *theColumns=(int*)omAlloc(k*sizeof(int));
577 mk.getAbsoluteColumnIndices(theColumns);
578 /* the next line provides the return value for the case k = 1 */
579 int e = getEntry(theRows[0], theColumns[0]);
580 if (characteristic != 0) e = e % characteristic;
581 if (iSB != 0) e = getReduction(e, iSB);
582 IntMinorValue mv(e, 0, 0, 0, 0, -1, -1);
583 if (k > 1)
584 {
585 /* the matrix to perform Bareiss with */
586 long *tempMatrix=(long*)omAlloc(k * k *sizeof(long));
587 /* copy correct set of entries from _intMatrix to tempMatrix */
588 int i = 0;
589 for (int r = 0; r < k; r++)
590 for (int c = 0; c < k; c++)
591 {
592 e = getEntry(theRows[r], theColumns[c]);
593 if (characteristic != 0) e = e % characteristic;
594 tempMatrix[i++] = e;
595 }
596 /* Bareiss algorithm operating on tempMatrix which is at least 2x2 */
597 int sign = 1; /* This will store the correct sign resulting
598 from permuting the rows of tempMatrix. */
599 int *rowPermutation=(int*)omAlloc(k*sizeof(int));
600 /* This is for storing the permutation of rows
601 resulting from searching for a non-zero
602 pivot element. */
603 for (int i = 0; i < k; i++) rowPermutation[i] = i;
604 int divisor = 1; /* the Bareiss divisor */
605 for (int r = 0; r <= k - 2; r++)
606 {
607 /* look for a non-zero entry in column r: */
608 int i = r;
609 while ((i < k) && (tempMatrix[rowPermutation[i] * k + r] == 0))
610 i++;
611 if (i == k)
612 /* There is no non-zero entry; hence the minor is zero. */
613 return IntMinorValue(0, 0, 0, 0, 0, -1, -1);
614 if (i != r)
615 {
616 /* We swap the rows with indices r and i: */
617 int j = rowPermutation[i];
618 rowPermutation[i] = rowPermutation[r];
619 rowPermutation[r] = j;
620 /* Now we know that tempMatrix[rowPermutation[r] * k + r] is not zero.
621 But careful; we have to negate the sign, as there is always an odd
622 number of row transpositions to swap two given rows of a matrix. */
623 sign = -sign;
624 }
625 if (r >= 1) divisor = tempMatrix[rowPermutation[r - 1] * k + r - 1];
626 for (int rr = r + 1; rr < k; rr++)
627 for (int cc = r + 1; cc < k; cc++)
628 {
629 e = rowPermutation[rr] * k + cc;
630 /* Attention: The following may cause an overflow and
631 thus a wrong result: */
632 tempMatrix[e] = tempMatrix[e] * tempMatrix[rowPermutation[r] * k + r]
633 - tempMatrix[rowPermutation[r] * k + cc]
634 * tempMatrix[rowPermutation[rr] * k + r];
635 /* The following is, by theory, always a division without
636 remainder: */
637 tempMatrix[e] = tempMatrix[e] / divisor;
638 if (characteristic != 0)
639 tempMatrix[e] = tempMatrix[e] % characteristic;
640 }
641 omFree(rowPermutation);
642 omFree(tempMatrix);
643 }
644 int theValue = sign * tempMatrix[rowPermutation[k - 1] * k + k - 1];
645 if (iSB != 0) theValue = getReduction(theValue, iSB);
646 mv = IntMinorValue(theValue, 0, 0, 0, 0, -1, -1);
647 }
648 omFree(theRows);
649 omFree(theColumns);
650 return mv;
651}
652
653int IntMinorProcessor::getEntry (const int rowIndex,
654 const int columnIndex) const
655{
656 return _intMatrix[rowIndex * _columns + columnIndex];
657}
658
660 const int k, const MinorKey& mk,
661 const bool multipleMinors,
663 const int characteristic, const ideal& iSB)
664{
665 assume(k > 0); /* k is the minor's dimension; the minor must be at least
666 1x1 */
667 /* The method works by recursion, and using Lapace's Theorem along
668 the row/column with the most zeros. */
669 if (k == 1)
670 {
672 if (characteristic != 0) e = e % characteristic;
673 if (iSB != 0) e = getReduction(e, iSB);
674 return IntMinorValue(e, 0, 0, 0, 0, -1, -1);
675 /* we set "-1" as, for k == 1, we do not have any cache retrievals */
676 }
677 else
678 {
679 int b = getBestLine(k, mk); /* row or column with
680 most zeros */
681 int result = 0; /* This will contain the
682 value of the minor. */
683 int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions
684 and multiplications,
685 ..."a*" for
686 accumulated operation
687 counters */
688 IntMinorValue mv(0, 0, 0, 0, 0, 0, 0); /* for storing all
689 intermediate minors */
690 bool hadNonZeroEntry = false;
691 if (b >= 0)
692 {
693 /* This means that the best line is the row with absolute (0-based)
694 index b.
695 Using Laplace, the sign of the contributing minors must be
696 iterating; the initial sign depends on the relative index of b
697 in minorRowKey: */
698 int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
699 for (int c = 0; c < k; c++) /* This iterates over all involved
700 columns. */
701 {
702 int absoluteC = mk.getAbsoluteColumnIndex(c);
703 if (getEntry(b, absoluteC) != 0) /* Only then do we have to consider
704 this sub-determinante. */
705 {
706 hadNonZeroEntry = true;
707 /* Next MinorKey is mk with row b and column absoluteC omitted. */
708 MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
709 if (cch.hasKey(subMk))
710 { /* trying to find the result in the cache */
711 mv = cch.getValue(subMk);
712 mv.incrementRetrievals(); /* once more, we made use of the cached
713 value for key mk */
714 cch.put(subMk, mv); /* We need to do this with "put", as the
715 (altered) number of retrievals may have
716 an impact on the internal ordering among
717 the cached entries. */
718 }
719 else
720 {
721 mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch,
722 characteristic, iSB); /* recursive call */
723 /* As this minor was not in the cache, we count the additions
724 and multiplications that we needed to perform in the
725 recursive call: */
726 m += mv.getMultiplications();
727 s += mv.getAdditions();
728 }
729 /* In any case, we count all nested operations in the accumulative
730 counters: */
732 as += mv.getAccumulatedAdditions();
733 /* adding sub-determinante times matrix entry times appropriate
734 sign */
735 result += sign * mv.getResult() * getEntry(b, absoluteC);
736 if (characteristic != 0) result = result % characteristic;
737 s++; m++; as++; am++; /* This is for the last addition and
738 multiplication. */
739 }
740 sign = - sign; /* alternating the sign */
741 }
742 }
743 else
744 {
745 b = - b - 1;
746 /* This means that the best line is the column with absolute (0-based)
747 index b.
748 Using Laplace, the sign of the contributing minors must be iterating;
749 the initial sign depends on the relative index of b in
750 minorColumnKey: */
751 int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
752 for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
753 {
754 int absoluteR = mk.getAbsoluteRowIndex(r);
755 if (getEntry(absoluteR, b) != 0) /* Only then do we have to consider
756 this sub-determinante. */
757 {
758 hadNonZeroEntry = true;
759 /* Next MinorKey is mk with row absoluteR and column b omitted. */
760 MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
761 if (cch.hasKey(subMk))
762 { /* trying to find the result in the cache */
763 mv = cch.getValue(subMk);
764 mv.incrementRetrievals(); /* once more, we made use of the cached
765 value for key mk */
766 cch.put(subMk, mv); /* We need to do this with "put", as the
767 (altered) number of retrievals may have an
768 impact on the internal ordering among the
769 cached entries. */
770 }
771 else
772 {
773 mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch, characteristic, iSB); /* recursive call */
774 /* As this minor was not in the cache, we count the additions and
775 multiplications that we needed to do in the recursive call: */
776 m += mv.getMultiplications();
777 s += mv.getAdditions();
778 }
779 /* In any case, we count all nested operations in the accumulative
780 counters: */
782 as += mv.getAccumulatedAdditions();
783 /* adding sub-determinante times matrix entry times appropriate
784 sign: */
785 result += sign * mv.getResult() * getEntry(absoluteR, b);
786 if (characteristic != 0) result = result % characteristic;
787 s++; m++; as++; am++; /* This is for the last addition and
788 multiplication. */
789 }
790 sign = - sign; /* alternating the sign */
791 }
792 }
793 /* Let's cache the newly computed minor: */
794 int potentialRetrievals = NumberOfRetrievals(_containerRows,
796 _minorSize, k,
797 multipleMinors);
798 if (hadNonZeroEntry)
799 {
800 s--; as--; /* first addition was 0 + ..., so we do not count it */
801 }
802 if (s < 0) s = 0; /* may happen when all subminors are zero and no
803 addition needs to be performed */
804 if (as < 0) as = 0; /* may happen when all subminors are zero and no
805 addition needs to be performed */
806 if (iSB != 0) result = getReduction(result, iSB);
807 IntMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
808 cch.put(mk, newMV); /* Here's the actual put inside the cache. */
809 return newMV;
810 }
811}
812
814{
816}
817
818poly PolyMinorProcessor::getEntry (const int rowIndex,
819 const int columnIndex) const
820{
821 return _polyMatrix[rowIndex * _columns + columnIndex];
822}
823
824bool PolyMinorProcessor::isEntryZero (const int absoluteRowIndex,
825 const int absoluteColumnIndex) const
826{
827 return getEntry(absoluteRowIndex, absoluteColumnIndex) == NULL;
828}
829
831{
832 char h[32];
833 string t = "";
834 string s = "PolyMinorProcessor:";
835 s += "\n matrix: ";
836 sprintf(h, "%d", _rows); s += h;
837 s += " x ";
838 sprintf(h, "%d", _columns); s += h;
839 int myIndexArray[500];
840 s += "\n considered submatrix has row indices: ";
841 _container.getAbsoluteRowIndices(myIndexArray);
842 for (int k = 0; k < _containerRows; k++)
843 {
844 if (k != 0) s += ", ";
845 sprintf(h, "%d", myIndexArray[k]); s += h;
846 }
847 s += " (first row of matrix has index 0)";
848 s += "\n considered submatrix has column indices: ";
850 for (int k = 0; k < _containerColumns; k++)
851 {
852 if (k != 0) s += ", ";
853 sprintf(h, "%d", myIndexArray[k]); s += h;
854 }
855 s += " (first column of matrix has index 0)";
856 s += "\n size of considered minor(s): ";
857 sprintf(h, "%d", _minorSize); s += h;
858 s += "x";
859 s += h;
860 return s;
861}
862
864{
865 /* free memory of _polyMatrix */
866 int n = _rows * _columns;
867 for (int i = 0; i < n; i++)
870}
871
872void PolyMinorProcessor::defineMatrix (const int numberOfRows,
873 const int numberOfColumns,
874 const poly* polyMatrix)
875{
876 /* free memory of _polyMatrix */
877 int n = _rows * _columns;
878 for (int i = 0; i < n; i++)
881
882 _rows = numberOfRows;
883 _columns = numberOfColumns;
884 n = _rows * _columns;
885
886 /* allocate memory for new entries in _polyMatrix */
887 _polyMatrix = (poly*)omAlloc(n*sizeof(poly));
888
889 /* copying values from one-dimensional method
890 parameter "polyMatrix" */
891 for (int i = 0; i < n; i++)
892 _polyMatrix[i] = pCopy(polyMatrix[i]);
893}
894
896 const int* rowIndices,
897 const int* columnIndices,
899 const ideal& iSB)
900{
901 defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
903 /* call a helper method which recursively computes the minor using the cache
904 c: */
905 return getMinorPrivateLaplace(dimension, _container, false, c, iSB);
906}
907
909 const int* rowIndices,
910 const int* columnIndices,
911 const char* algorithm,
912 const ideal& iSB)
913{
914 defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
916 /* call a helper method which computes the minor (without using a cache): */
917 if (strcmp(algorithm, "Laplace") == 0)
919 else if (strcmp(algorithm, "Bareiss") == 0)
921 else assume(false);
922
923 /* The following code is never reached and just there to make the
924 compiler happy: */
925 return PolyMinorValue();
926}
927
929 const ideal& iSB)
930{
931 /* call a helper method which computes the minor (without using a
932 cache): */
933 if (strcmp(algorithm, "Laplace") == 0)
935 else if (strcmp(algorithm, "Bareiss") == 0)
937 else assume(false);
938
939 /* The following code is never reached and just there to make the
940 compiler happy: */
941 return PolyMinorValue();
942}
943
945 PolyMinorValue>& c,
946 const ideal& iSB)
947{
948 /* computation with cache */
949 return getMinorPrivateLaplace(_minorSize, _minor, true, c, iSB);
950}
951
952/* assumes that all entries in polyMatrix are already reduced w.r.t. iSB */
954 const MinorKey& mk,
955 const ideal& iSB)
956{
957 assume(k > 0); /* k is the minor's dimension; the minor must be at least
958 1x1 */
959 /* The method works by recursion, and using Lapace's Theorem along the
960 row/column with the most zeros. */
961 if (k == 1)
962 {
965 0, 0, 0, 0, -1, -1);
966 /* "-1" is to signal that any statistics about the number of retrievals
967 does not make sense, as we do not use a cache. */
968 return pmv;
969 }
970 else
971 {
972 /* Here, the minor must be 2x2 or larger. */
973 int b = getBestLine(k, mk); /* row or column with most
974 zeros */
975 poly result = NULL; /* This will contain the
976 value of the minor. */
977 int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions
978 and multiplications,
979 ..."a*" for accumulated
980 operation counters */
981 bool hadNonZeroEntry = false;
982 if (b >= 0)
983 {
984 /* This means that the best line is the row with absolute (0-based)
985 index b.
986 Using Laplace, the sign of the contributing minors must be iterating;
987 the initial sign depends on the relative index of b in minorRowKey: */
988 int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
989 for (int c = 0; c < k; c++) /* This iterates over all involved columns. */
990 {
991 int absoluteC = mk.getAbsoluteColumnIndex(c);
992 if (!isEntryZero(b, absoluteC)) /* Only then do we have to consider
993 this sub-determinante. */
994 {
995 hadNonZeroEntry = true;
996 /* Next MinorKey is mk with row b and column absoluteC omitted. */
997 MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
998 /* recursive call: */
999 PolyMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, iSB);
1000 m += mv.getMultiplications();
1001 s += mv.getAdditions();
1003 as += mv.getAccumulatedAdditions();
1004 poly signPoly = pISet(sign);
1005 poly temp = pp_Mult_qq(mv.getResult(), getEntry(b, absoluteC),
1006 currRing);
1007 temp = p_Mult_q(signPoly, temp, currRing);
1008 result = p_Add_q(result, temp, currRing);
1009#ifdef COUNT_AND_PRINT_OPERATIONS
1010 multsPoly++;
1011 addsPoly++;
1012 multsMon += pLength(mv.getResult()) * pLength(getEntry(b, absoluteC));
1013#endif
1014 s++; m++; as++, am++; /* This is for the addition and multiplication
1015 in the previous lines of code. */
1016 }
1017 sign = - sign; /* alternating the sign */
1018 }
1019 }
1020 else
1021 {
1022 b = - b - 1;
1023 /* This means that the best line is the column with absolute (0-based)
1024 index b.
1025 Using Laplace, the sign of the contributing minors must be iterating;
1026 the initial sign depends on the relative index of b in
1027 minorColumnKey: */
1028 int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
1029 for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
1030 {
1031 int absoluteR = mk.getAbsoluteRowIndex(r);
1032 if (!isEntryZero(absoluteR, b)) /* Only then do we have to consider
1033 this sub-determinante. */
1034 {
1035 hadNonZeroEntry = true;
1036 /* This is mk with row absoluteR and column b omitted. */
1037 MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
1038 /* recursive call: */
1039 PolyMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, iSB);
1040 m += mv.getMultiplications();
1041 s += mv.getAdditions();
1043 as += mv.getAccumulatedAdditions();
1044 poly signPoly = pISet(sign);
1045 poly temp = pp_Mult_qq(mv.getResult(), getEntry(absoluteR, b),
1046 currRing);
1047 temp = p_Mult_q(signPoly, temp, currRing);
1048 result = p_Add_q(result, temp, currRing);
1049#ifdef COUNT_AND_PRINT_OPERATIONS
1050 multsPoly++;
1051 addsPoly++;
1052 multsMon += pLength(mv.getResult()) * pLength(getEntry(absoluteR, b));
1053#endif
1054 s++; m++; as++, am++; /* This is for the addition and multiplication
1055 in the previous lines of code. */
1056 }
1057 sign = - sign; /* alternating the sign */
1058 }
1059 }
1060 if (hadNonZeroEntry)
1061 {
1062 s--; as--; /* first addition was 0 + ..., so we do not count it */
1063#ifdef COUNT_AND_PRINT_OPERATIONS
1064 addsPoly--;
1065#endif
1066 }
1067 if (s < 0) s = 0; /* may happen when all subminors are zero and no
1068 addition needs to be performed */
1069 if (as < 0) as = 0; /* may happen when all subminors are zero and no
1070 addition needs to be performed */
1071 if (iSB != NULL)
1072 {
1073 poly tmpresult = kNF(iSB, currRing->qideal, result);
1074 pDelete(&result);
1075 result=tmpresult;
1076 }
1077 PolyMinorValue newMV(result, m, s, am, as, -1, -1);
1078 /* "-1" is to signal that any statistics about the number of retrievals
1079 does not make sense, as we do not use a cache. */
1080 pDelete(&result);
1081 return newMV;
1082 }
1083}
1084
1085/* assumes that all entries in polyMatrix are already reduced w.r.t. iSB */
1087 const int k,
1088 const MinorKey& mk,
1089 const bool multipleMinors,
1091 const ideal& iSB)
1092{
1093 assume(k > 0); /* k is the minor's dimension; the minor must be at least
1094 1x1 */
1095 /* The method works by recursion, and using Lapace's Theorem along
1096 the row/column with the most zeros. */
1097 if (k == 1)
1098 {
1101 0, 0, 0, 0, -1, -1);
1102 /* we set "-1" as, for k == 1, we do not have any cache retrievals */
1103 return pmv;
1104 }
1105 else
1106 {
1107 int b = getBestLine(k, mk); /* row or column with most
1108 zeros */
1109 poly result = NULL; /* This will contain the
1110 value of the minor. */
1111 int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions
1112 and multiplications,
1113 ..."a*" for accumulated
1114 operation counters */
1115 bool hadNonZeroEntry = false;
1116 if (b >= 0)
1117 {
1118 /* This means that the best line is the row with absolute (0-based)
1119 index b.
1120 Using Laplace, the sign of the contributing minors must be iterating;
1121 the initial sign depends on the relative index of b in
1122 minorRowKey: */
1123 int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
1124 for (int c = 0; c < k; c++) /* This iterates over all involved columns. */
1125 {
1126 int absoluteC = mk.getAbsoluteColumnIndex(c);
1127 if (!isEntryZero(b, absoluteC)) /* Only then do we have to consider
1128 this sub-determinante. */
1129 {
1130 hadNonZeroEntry = true;
1131 PolyMinorValue mv; /* for storing all intermediate minors */
1132 /* Next MinorKey is mk with row b and column absoluteC omitted. */
1133 MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
1134 if (cch.hasKey(subMk))
1135 { /* trying to find the result in the cache */
1136 mv = cch.getValue(subMk);
1137 mv.incrementRetrievals(); /* once more, we made use of the cached
1138 value for key mk */
1139 cch.put(subMk, mv); /* We need to do this with "put", as the
1140 (altered) number of retrievals may have an
1141 impact on the internal ordering among cache
1142 entries. */
1143 }
1144 else
1145 {
1146 /* recursive call: */
1147 mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch,
1148 iSB);
1149 /* As this minor was not in the cache, we count the additions and
1150 multiplications that we needed to do in the recursive call: */
1151 m += mv.getMultiplications();
1152 s += mv.getAdditions();
1153 }
1154 /* In any case, we count all nested operations in the accumulative
1155 counters: */
1157 as += mv.getAccumulatedAdditions();
1158 poly signPoly = pISet(sign);
1159 poly temp = pp_Mult_qq(mv.getResult(), getEntry(b, absoluteC),
1160 currRing);
1161 temp = p_Mult_q(signPoly, temp, currRing);
1162 result = p_Add_q(result, temp, currRing);
1163#ifdef COUNT_AND_PRINT_OPERATIONS
1164 multsPoly++;
1165 addsPoly++;
1166 multsMon += pLength(mv.getResult()) * pLength(getEntry(b, absoluteC));
1167#endif
1168 s++; m++; as++; am++; /* This is for the addition and multiplication
1169 in the previous lines of code. */
1170 }
1171 sign = - sign; /* alternating the sign */
1172 }
1173 }
1174 else
1175 {
1176 b = - b - 1;
1177 /* This means that the best line is the column with absolute (0-based)
1178 index b.
1179 Using Laplace, the sign of the contributing minors must be iterating;
1180 the initial sign depends on the relative index of b in
1181 minorColumnKey: */
1182 int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
1183 for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
1184 {
1185 int absoluteR = mk.getAbsoluteRowIndex(r);
1186 if (!isEntryZero(absoluteR, b)) /* Only then do we have to consider
1187 this sub-determinante. */
1188 {
1189 hadNonZeroEntry = true;
1190 PolyMinorValue mv; /* for storing all intermediate minors */
1191 /* Next MinorKey is mk with row absoluteR and column b omitted. */
1192 MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
1193 if (cch.hasKey(subMk))
1194 { /* trying to find the result in the cache */
1195 mv = cch.getValue(subMk);
1196 mv.incrementRetrievals(); /* once more, we made use of the cached
1197 value for key mk */
1198 cch.put(subMk, mv); /* We need to do this with "put", as the
1199 (altered) number of retrievals may have an
1200 impact on the internal ordering among the
1201 cached entries. */
1202 }
1203 else
1204 {
1205 mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch,
1206 iSB); /* recursive call */
1207 /* As this minor was not in the cache, we count the additions and
1208 multiplications that we needed to do in the recursive call: */
1209 m += mv.getMultiplications();
1210 s += mv.getAdditions();
1211 }
1212 /* In any case, we count all nested operations in the accumulative
1213 counters: */
1215 as += mv.getAccumulatedAdditions();
1216 poly signPoly = pISet(sign);
1217 poly temp = pp_Mult_qq(mv.getResult(), getEntry(absoluteR, b),
1218 currRing);
1219 temp = p_Mult_q(signPoly, temp, currRing);
1220 result = p_Add_q(result, temp, currRing);
1221#ifdef COUNT_AND_PRINT_OPERATIONS
1222 multsPoly++;
1223 addsPoly++;
1224 multsMon += pLength(mv.getResult()) * pLength(getEntry(absoluteR, b));
1225#endif
1226 s++; m++; as++; am++; /* This is for the addition and multiplication
1227 in the previous lines of code. */
1228 }
1229 sign = - sign; /* alternating the sign */
1230 }
1231 }
1232 /* Let's cache the newly computed minor: */
1233 int potentialRetrievals = NumberOfRetrievals(_containerRows,
1235 _minorSize,
1236 k,
1237 multipleMinors);
1238 if (hadNonZeroEntry)
1239 {
1240 s--; as--; /* first addition was 0 + ..., so we do not count it */
1241#ifdef COUNT_AND_PRINT_OPERATIONS
1242 addsPoly--;
1243#endif
1244 }
1245 if (s < 0) s = 0; /* may happen when all subminors are zero and no
1246 addition needs to be performed */
1247 if (as < 0) as = 0; /* may happen when all subminors are zero and no
1248 addition needs to be performed */
1249 if (iSB != NULL)
1250 {
1251 poly tmpresult = kNF(iSB, currRing->qideal, result);
1252 pDelete(&tmpresult);
1253 result=tmpresult;
1254 }
1255 PolyMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
1257 cch.put(mk, newMV); /* Here's the actual put inside the cache. */
1258 return newMV;
1259 }
1260}
1261
1262/* This can only be used in the case of coefficients coming from a field
1263 or at least an integral domain. */
1264static void addOperationBucket(poly f1, poly f2, kBucket_pt bucket)
1265{
1266 /* fills all terms of f1 * f2 into the bucket */
1267 poly a = f1; poly b = f2;
1268 int aLen = pLength(a); int bLen = pLength(b);
1269 if (aLen > bLen)
1270 {
1271 b = f1; a = f2; bLen = aLen;
1272 }
1273 pNormalize(b);
1274
1275 while (a != NULL)
1276 {
1277 /* The next line actually uses only LT(a): */
1278 kBucket_Plus_mm_Mult_pp(bucket, a, b, bLen);
1279 a = pNext(a);
1280 }
1281}
1282
1283/* computes the polynomial (p1 * p2 - p3 * p4) and puts result into p1;
1284 the method destroys the old value of p1;
1285 p2, p3, and p4 may be pNormalize-d but must, apart from that,
1286 not be changed;
1287 This can only be used in the case of coefficients coming from a field
1288 or at least an integral domain. */
1289static void elimOperationBucketNoDiv(poly &p1, poly p2, poly p3, poly p4)
1290{
1291#ifdef COUNT_AND_PRINT_OPERATIONS
1292 if ((pLength(p1) != 0) && (pLength(p2) != 0))
1293 {
1294 multsPoly++;
1295 multsMon += pLength(p1) * pLength(p2);
1296 }
1297 if ((pLength(p3) != 0) && (pLength(p4) != 0))
1298 {
1299 multsPoly++;
1300 multsMon += pLength(p3) * pLength(p4);
1301 }
1302 if ((pLength(p1) != 0) && (pLength(p2) != 0) &&
1303 (pLength(p3) != 0) && (pLength(p4) != 0))
1304 addsPoly++;
1305#endif
1306 kBucket_pt myBucket = kBucketCreate(currRing);
1307 addOperationBucket(p1, p2, myBucket);
1308 poly p3Neg = pNeg(pCopy(p3));
1309 addOperationBucket(p3Neg, p4, myBucket);
1310 pDelete(&p3Neg);
1311 pDelete(&p1);
1312 p1 = kBucketClear(myBucket);
1313 kBucketDestroy(&myBucket);
1314}
1315
1316/* computes the polynomial (p1 * p2 - p3 * p4) / p5 and puts result into p1;
1317 the method destroys the old value of p1;
1318 p2, p3, p4, and p5 may be pNormalize-d but must, apart from that,
1319 not be changed;
1320 c5 is assumed to be the leading coefficient of p5;
1321 p5Len is assumed to be the length of p5;
1322 This can only be used in the case of coefficients coming from a field
1323 or at least an integral domain. */
1324void elimOperationBucket(poly &p1, poly &p2, poly &p3, poly &p4, poly &p5,
1325 number &c5, int p5Len)
1326{
1327#ifdef COUNT_AND_PRINT_OPERATIONS
1328 if ((pLength(p1) != 0) && (pLength(p2) != 0))
1329 {
1330 multsPoly++;
1331 multsMon += pLength(p1) * pLength(p2);
1332 }
1333 if ((pLength(p3) != 0) && (pLength(p4) != 0))
1334 {
1335 multsPoly++;
1336 multsMon += pLength(p3) * pLength(p4);
1337 }
1338 if ((pLength(p1) != 0) && (pLength(p2) != 0) &&
1339 (pLength(p3) != 0) && (pLength(p4) != 0))
1340 addsPoly++;
1341#endif
1342 kBucket_pt myBucket = kBucketCreate(currRing);
1343 addOperationBucket(p1, p2, myBucket);
1344 poly p3Neg = pNeg(pCopy(p3));
1345 addOperationBucket(p3Neg, p4, myBucket);
1346 pDelete(&p3Neg);
1347
1348 /* Now, myBucket contains all terms of p1 * p2 - p3 * p4.
1349 Now we need to perform the polynomial division myBucket / p5
1350 which is known to work without remainder: */
1351 pDelete(&p1); poly helperPoly = NULL;
1352
1353 poly bucketLm = pCopy(kBucketGetLm(myBucket));
1354 while (bucketLm != NULL)
1355 {
1356 /* divide bucketLm by the leading term of p5 and put result into bucketLm;
1357 we start with the coefficients;
1358 note that bucketLm will always represent a term */
1359 number coeff = nDiv(pGetCoeff(bucketLm), c5);
1360 nNormalize(coeff);
1361 pSetCoeff(bucketLm, coeff);
1362 /* subtract exponent vector of p5 from that of quotient; modifies
1363 quotient */
1364 p_ExpVectorSub(bucketLm, p5, currRing);
1365#ifdef COUNT_AND_PRINT_OPERATIONS
1366 divsMon++;
1367 multsMonForDiv += p5Len;
1368 multsMon += p5Len;
1369 savedMultsMFD++;
1370 multsPoly++;
1371 multsPolyForDiv++;
1372 addsPoly++;
1373 addsPolyForDiv++;
1374#endif
1375 kBucket_Minus_m_Mult_p(myBucket, bucketLm, p5, &p5Len);
1376 /* The following lines make bucketLm the new leading term of p1,
1377 i.e., put bucketLm in front of everything which is already in p1.
1378 Thus, after the while loop, we need to revert p1. */
1379 helperPoly = bucketLm;
1380 helperPoly->next = p1;
1381 p1 = helperPoly;
1382
1383 bucketLm = pCopy(kBucketGetLm(myBucket));
1384 }
1385 p1 = pReverse(p1);
1386 kBucketDestroy(&myBucket);
1387}
1388
1389/* assumes that all entries in polyMatrix are already reduced w.r.t. iSB
1390 This can only be used in the case of coefficients coming from a field!!! */
1392 const MinorKey& mk,
1393 const ideal& iSB)
1394{
1395 assume(k > 0); /* k is the minor's dimension; the minor must be at least
1396 1x1 */
1397 int *theRows=(int*)omAlloc(k*sizeof(int));
1398 mk.getAbsoluteRowIndices(theRows);
1399 int *theColumns=(int*)omAlloc(k*sizeof(int));
1400 mk.getAbsoluteColumnIndices(theColumns);
1401 if (k == 1)
1402 {
1403 PolyMinorValue tmp=PolyMinorValue(getEntry(theRows[0], theColumns[0]),
1404 0, 0, 0, 0, -1, -1);
1405 omFree(theColumns);
1406 omFree(theRows);
1407 return tmp;
1408 }
1409 else /* k > 0 */
1410 {
1411 /* the matrix to perform Bareiss with */
1412 poly* tempMatrix = (poly*)omAlloc(k * k * sizeof(poly));
1413 /* copy correct set of entries from _polyMatrix to tempMatrix */
1414 int i = 0;
1415 for (int r = 0; r < k; r++)
1416 for (int c = 0; c < k; c++)
1417 tempMatrix[i++] = pCopy(getEntry(theRows[r], theColumns[c]));
1418
1419 /* Bareiss algorithm operating on tempMatrix which is at least 2x2 */
1420 int sign = 1; /* This will store the correct sign resulting from
1421 permuting the rows of tempMatrix. */
1422 int *rowPermutation=(int*)omAlloc(k*sizeof(int));
1423 /* This is for storing the permutation of rows
1424 resulting from searching for a non-zero pivot
1425 element. */
1426 for (int i = 0; i < k; i++) rowPermutation[i] = i;
1427 poly divisor = NULL;
1428 int divisorLength = 0;
1429 number divisorLC;
1430 for (int r = 0; r <= k - 2; r++)
1431 {
1432 /* look for a non-zero entry in column r, rows = r .. (k - 1)
1433 s.t. the polynomial has least complexity: */
1434 int minComplexity = -1; int complexity = 0; int bestRow = -1;
1435 poly pp = NULL;
1436 for (int i = r; i < k; i++)
1437 {
1438 pp = tempMatrix[rowPermutation[i] * k + r];
1439 if (pp != NULL)
1440 {
1441 if (minComplexity == -1)
1442 {
1443 minComplexity = pSize(pp);
1444 bestRow = i;
1445 }
1446 else
1447 {
1448 complexity = 0;
1449 while ((pp != NULL) && (complexity < minComplexity))
1450 {
1451 complexity += nSize(pGetCoeff(pp)); pp = pNext(pp);
1452 }
1453 if (complexity < minComplexity)
1454 {
1455 minComplexity = complexity;
1456 bestRow = i;
1457 }
1458 }
1459 if (minComplexity <= 1) break; /* terminate the search */
1460 }
1461 }
1462 if (bestRow == -1)
1463 {
1464 /* There is no non-zero entry; hence the minor is zero. */
1465 for (int i = 0; i < k * k; i++) pDelete(&tempMatrix[i]);
1466 return PolyMinorValue(NULL, 0, 0, 0, 0, -1, -1);
1467 }
1468 pNormalize(tempMatrix[rowPermutation[bestRow] * k + r]);
1469 if (bestRow != r)
1470 {
1471 /* We swap the rows with indices r and i: */
1472 int j = rowPermutation[bestRow];
1473 rowPermutation[bestRow] = rowPermutation[r];
1474 rowPermutation[r] = j;
1475 /* Now we know that tempMatrix[rowPermutation[r] * k + r] is not zero.
1476 But careful; we have to negate the sign, as there is always an odd
1477 number of row transpositions to swap two given rows of a matrix. */
1478 sign = -sign;
1479 }
1480#if (defined COUNT_AND_PRINT_OPERATIONS) && (COUNT_AND_PRINT_OPERATIONS > 2)
1481 poly w = NULL; int wl = 0;
1482 printf("matrix after %d steps:\n", r);
1483 for (int u = 0; u < k; u++)
1484 {
1485 for (int v = 0; v < k; v++)
1486 {
1487 if ((v < r) && (u > v))
1488 wl = 0;
1489 else
1490 {
1491 w = tempMatrix[rowPermutation[u] * k + v];
1492 wl = pLength(w);
1493 }
1494 printf("%5d ", wl);
1495 }
1496 printf("\n");
1497 }
1498 printCounters ("", false);
1499#endif
1500 if (r != 0)
1501 {
1502 divisor = tempMatrix[rowPermutation[r - 1] * k + r - 1];
1503 pNormalize(divisor);
1504 divisorLength = pLength(divisor);
1505 divisorLC = pGetCoeff(divisor);
1506 }
1507 for (int rr = r + 1; rr < k; rr++)
1508 for (int cc = r + 1; cc < k; cc++)
1509 {
1510 if (r == 0)
1511 elimOperationBucketNoDiv(tempMatrix[rowPermutation[rr] * k + cc],
1512 tempMatrix[rowPermutation[r] * k + r],
1513 tempMatrix[rowPermutation[r] * k + cc],
1514 tempMatrix[rowPermutation[rr] * k + r]);
1515 else
1516 elimOperationBucket(tempMatrix[rowPermutation[rr] * k + cc],
1517 tempMatrix[rowPermutation[r] * k + r],
1518 tempMatrix[rowPermutation[r] * k + cc],
1519 tempMatrix[rowPermutation[rr] * k + r],
1520 divisor, divisorLC, divisorLength);
1521 }
1522 }
1523#if (defined COUNT_AND_PRINT_OPERATIONS) && (COUNT_AND_PRINT_OPERATIONS > 2)
1524 poly w = NULL; int wl = 0;
1525 printf("matrix after %d steps:\n", k - 1);
1526 for (int u = 0; u < k; u++)
1527 {
1528 for (int v = 0; v < k; v++)
1529 {
1530 if ((v < k - 1) && (u > v))
1531 wl = 0;
1532 else
1533 {
1534 w = tempMatrix[rowPermutation[u] * k + v];
1535 wl = pLength(w);
1536 }
1537 printf("%5d ", wl);
1538 }
1539 printf("\n");
1540 }
1541#endif
1542 poly result = tempMatrix[rowPermutation[k - 1] * k + k - 1];
1543 tempMatrix[rowPermutation[k - 1] * k + k - 1]=NULL; /*value now in result*/
1544 if (sign == -1) result = pNeg(result);
1545 if (iSB != NULL)
1546 {
1547 poly tmpresult = kNF(iSB, currRing->qideal, result);
1548 pDelete(&result);
1549 result=tmpresult;
1550 }
1551 PolyMinorValue mv(result, 0, 0, 0, 0, -1, -1);
1552 for (int i = 0; i < k * k; i++) pDelete(&tempMatrix[i]);
1553 omFreeSize(tempMatrix, k * k * sizeof(poly));
1554 omFree(rowPermutation);
1555 omFree(theColumns);
1556 omFree(theRows);
1557 return mv;
1558 }
1559}
1560
static void addOperationBucket(poly f1, poly f2, kBucket_pt bucket)
int getReduction(const int i, const ideal &iSB)
static void elimOperationBucketNoDiv(poly &p1, poly p2, poly p3, poly p4)
void elimOperationBucket(poly &p1, poly &p2, poly &p3, poly &p4, poly &p5, number &c5, int p5Len)
void printCounters(char *prefix, bool resetToZero)
BOOLEAN dimension(leftv res, leftv args)
Definition: bbcone.cc:757
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
FILE * f
Definition: checklibs.c:9
Class Cache is a template-implementation of a cache with arbitrary classes for representing keys and ...
Definition: Cache.h:69
bool put(const KeyClass &key, const ValueClass &value)
Inserts the pair (key --> value) in the cache.
bool hasKey(const KeyClass &key) const
Checks whether the cache contains a pair (k --> v) such that k equals the given key.
ValueClass getValue(const KeyClass &key) const
Returns the value for a given key.
std::string toString() const
A method for providing a printable version of the represented MinorProcessor.
~IntMinorProcessor()
A destructor for deleting an instance.
IntMinorProcessor()
A constructor for creating an instance.
int * _intMatrix
private store for integer matrix entries
IntMinorValue getMinorPrivateBareiss(const int k, const MinorKey &mk, const int characteristic, const ideal &iSB)
A method for computing the value of a minor using Bareiss's algorithm.
bool isEntryZero(const int absoluteRowIndex, const int absoluteColumnIndex) const
A method for testing whether a matrix entry is zero.
void defineMatrix(const int numberOfRows, const int numberOfColumns, const int *matrix)
A method for defining a matrix with integer entries.
IntMinorValue getMinor(const int dimension, const int *rowIndices, const int *columnIndices, const int characteristic, const ideal &iSB, const char *algorithm)
A method for computing the value of a minor without using a cache.
IntMinorValue getNextMinor(const int characteristic, const ideal &iSB, const char *algorithm)
A method for obtaining the next minor when iterating through all minors of a given size within a pre-...
IntMinorValue getMinorPrivateLaplace(const int k, const MinorKey &mk, const bool multipleMinors, Cache< MinorKey, IntMinorValue > &c, int characteristic, const ideal &iSB)
A method for computing the value of a minor, using a cache.
int getEntry(const int rowIndex, const int columnIndex) const
A method for retrieving the matrix entry.
Class IntMinorValue is derived from MinorValue and can be used for representing values in a cache for...
Definition: Minor.h:718
int getResult() const
Accessor for the private field _result.
Definition: Minor.cc:1019
Class MinorKey can be used for representing keys in a cache for sub-determinantes; see class Cache.
Definition: Minor.h:40
void getAbsoluteColumnIndices(int *const target) const
A method for retrieving the 0-based indices of all columns encoded in this MinorKey.
Definition: Minor.cc:202
void selectFirstRows(const int k, const MinorKey &mk)
This method redefines the set of rows represented by this MinorKey.
Definition: Minor.cc:457
bool selectNextColumns(const int k, const MinorKey &mk)
This method redefines the set of columns represented by this MinorKey.
Definition: Minor.cc:669
void reset()
A method for deleting all entries of _rowKey and _columnKey.
Definition: Minor.cc:13
void selectFirstColumns(const int k, const MinorKey &mk)
This method redefines the set of columns represented by this MinorKey.
Definition: Minor.cc:498
MinorKey getSubMinorKey(const int absoluteEraseRowIndex, const int absoluteEraseColumnIndex) const
A method for retrieving a sub-MinorKey resulting from omitting one row and one column of this MinorKe...
Definition: Minor.cc:343
int getAbsoluteColumnIndex(const int i) const
A method for retrieving the (0-based) index of the i-th column in the set of columns encoded in this.
Definition: Minor.cc:149
void set(const int lengthOfRowArray, const unsigned int *rowKey, const int lengthOfColumnArray, const unsigned int *columnKey)
A setter method for class MinorKey.
Definition: Minor.cc:62
void getAbsoluteRowIndices(int *const target) const
A method for retrieving the 0-based indices of all rows encoded in this MinorKey.
Definition: Minor.cc:181
int getRelativeRowIndex(const int i) const
A method for retrieving the (0-based) relative index of the i-th row in this MinorKey.
Definition: Minor.cc:223
bool selectNextRows(const int k, const MinorKey &mk)
This method redefines the set of rows represented by this MinorKey.
Definition: Minor.cc:538
int getRelativeColumnIndex(const int i) const
A method for retrieving the (0-based) relative index of the i-th column in this MinorKey.
Definition: Minor.cc:255
int compare(const MinorKey &mk) const
A comparator for two instances of MinorKey.
Definition: Minor.cc:412
int getAbsoluteRowIndex(const int i) const
A method for retrieving the (0-based) index of the i-th row in the set of rows encoded in this.
Definition: Minor.cc:117
virtual bool isEntryZero(const int absoluteRowIndex, const int absoluteColumnIndex) const
A method for testing whether a matrix entry is zero.
static int IOverJ(const int i, const int j)
A static method for computing the binomial coefficient i over j.
MinorKey _minor
private store for the rows and columns of the minor of interest; Usually, this minor will encode subs...
void getCurrentColumnIndices(int *const target) const
A method for obtaining the current set of columns corresponding to the current minor when iterating t...
static int NumberOfRetrievals(const int rows, const int columns, const int containerMinorSize, const int minorSize, const bool multipleMinors)
A static method for computing the maximum number of retrievals of a minor.
static int Faculty(const int i)
A static method for computing the factorial of i.
void setMinorSize(const int minorSize)
Sets the size of the minor(s) of interest.
int _containerRows
private store for the number of rows in the container minor; This is set by MinorProcessor::defineSub...
virtual std::string toString() const
A method for providing a printable version of the represented MinorProcessor.
int getBestLine(const int k, const MinorKey &mk) const
A method for identifying the row or column with the most zeros.
bool setNextKeys(const int k)
A method for iterating through all possible subsets of k rows and k columns inside a pre-defined subm...
void print() const
A method for printing a string representation of the given MinorProcessor to std::cout.
int _columns
private store for the number of columns in the underlying matrix
int _minorSize
private store for the dimension of the minor(s) of interest
void defineSubMatrix(const int numberOfRows, const int *rowIndices, const int numberOfColumns, const int *columnIndices)
A method for defining a sub-matrix within a pre-defined matrix.
MinorKey _container
private store for the rows and columns of the container minor within the underlying matrix; _containe...
bool hasNextMinor()
A method for checking whether there is a next choice of rows and columns when iterating through all m...
virtual ~MinorProcessor()
A destructor for deleting an instance.
void getCurrentRowIndices(int *const target) const
A method for obtaining the current set of rows corresponding to the current minor when iterating thro...
MinorProcessor()
The default constructor.
int _rows
private store for the number of rows in the underlying matrix
int _containerColumns
private store for the number of columns in the container minor; This is set by MinorProcessor::define...
int getAdditions() const
A method for accessing the additions performed while computing this minor.
Definition: Minor.cc:888
int getAccumulatedAdditions() const
A method for accessing the additions performed while computing this minor, including all nested addit...
Definition: Minor.cc:898
int getMultiplications() const
A method for accessing the multiplications performed while computing this minor.
Definition: Minor.cc:883
void incrementRetrievals()
A method for incrementing the number of performed retrievals of this instance of MinorValue.
Definition: Minor.cc:873
int getAccumulatedMultiplications() const
A method for accessing the multiplications performed while computing this minor, including all nested...
Definition: Minor.cc:893
~PolyMinorProcessor()
A destructor for deleting an instance.
std::string toString() const
A method for providing a printable version of the represented MinorProcessor.
PolyMinorProcessor()
A constructor for creating an instance.
PolyMinorValue getNextMinor(const char *algorithm, const ideal &iSB)
A method for obtaining the next minor when iterating through all minors of a given size within a pre-...
bool isEntryZero(const int absoluteRowIndex, const int absoluteColumnIndex) const
A method for testing whether a matrix entry is zero.
poly getEntry(const int rowIndex, const int columnIndex) const
A method for retrieving the matrix entry.
void defineMatrix(const int numberOfRows, const int numberOfColumns, const poly *polyMatrix)
A method for defining a matrix with polynomial entries.
PolyMinorValue getMinor(const int dimension, const int *rowIndices, const int *columnIndices, const char *algorithm, const ideal &iSB)
A method for computing the value of a minor, without using a cache.
PolyMinorValue getMinorPrivateLaplace(const int k, const MinorKey &mk, const bool multipleMinors, Cache< MinorKey, PolyMinorValue > &c, const ideal &iSB)
A method for computing the value of a minor, using a cache.
PolyMinorValue getMinorPrivateBareiss(const int k, const MinorKey &mk, const ideal &iSB)
A method for computing the value of a minor, without using a cache.
poly * _polyMatrix
private store for polynomial matrix entries
Class PolyMinorValue is derived from MinorValue and can be used for representing values in a cache fo...
Definition: Minor.h:800
poly getResult() const
Accessor for the private field _result.
Definition: Minor.cc:1102
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition: coeffs.h:544
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm & w
Definition: facAbsFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
#define VAR
Definition: globaldefs.h:5
int binom(int n, int r)
STATIC_VAR int offset
Definition: janet.cc:29
STATIC_VAR Poly * h
Definition: janet.cc:971
void kBucketClear(kBucket_pt bucket, poly *p, int *length)
Definition: kbuckets.cc:521
void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l, poly spNoether)
Bpoly == Bpoly - m*p; where m is a monom Does not destroy p and m assume (*l <= 0 || pLength(p) == *l...
Definition: kbuckets.cc:722
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:216
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:209
void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
Bpoly == Bpoly + m*p; where m is a monom Does not destroy p and m assume (l <= 0 || pLength(p) == l)
Definition: kbuckets.cc:815
const poly kBucketGetLm(kBucket_pt bucket)
Definition: kbuckets.cc:506
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3182
#define assume(x)
Definition: mod2.h:389
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nDiv(a, b)
Definition: numbers.h:32
#define nSize(n)
Definition: numbers.h:39
#define nNormalize(n)
Definition: numbers.h:30
#define omfree(addr)
Definition: omAllocDecl.h:237
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
static int pLength(poly a)
Definition: p_polys.h:188
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1112
static void p_ExpVectorSub(poly p1, poly p2, const ring r)
Definition: p_polys.h:1438
static poly pReverse(poly p)
Definition: p_polys.h:333
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1149
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatibility layer for legacy polynomial operations (over currRing)
#define pDelete(p_ptr)
Definition: polys.h:186
#define pNeg(p)
Definition: polys.h:198
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pNormalize(p)
Definition: polys.h:317
#define pSize(p)
Definition: polys.h:318
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pISet(i)
Definition: polys.h:312
void PrintS(const char *s)
Definition: reporter.cc:284
static int sign(int x)
Definition: ring.cc:3427