My Project
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"
12 #include "kernel/GBEngine/kstd1.h"
13 
14 #include "kernel/ideals.h"
15 
16 using namespace std;
17 
18 #ifdef COUNT_AND_PRINT_OPERATIONS
19 VAR long addsPoly = 0; /* for the number of additions of two polynomials */
20 VAR long multsPoly = 0; /* for the number of multiplications of two polynomials */
21 VAR long addsPolyForDiv = 0; /* for the number of additions of two polynomials for
22  polynomial division part */
23 VAR long multsPolyForDiv = 0; /* for the number of multiplications of two polynomials
24  for polynomial division part */
25 VAR long multsMon = 0; /* for the number of multiplications of two monomials */
26 VAR long multsMonForDiv = 0; /* for the number of m-m-multiplications for polynomial
27  division part */
28 VAR 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) */
33 VAR 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) */
36 void 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 
57 int 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 
108 void MinorProcessor::setMinorSize(const int minorSize)
109 {
110  _minorSize = minorSize;
111  _minor.reset();
112 }
113 
115 {
116  return setNextKeys(_minorSize);
117 }
118 
119 void MinorProcessor::getCurrentRowIndices(int* const target) const
120 {
121  return _minor.getAbsoluteRowIndices(target);
122 }
123 
124 void MinorProcessor::getCurrentColumnIndices(int* const target) const
125 {
126  return _minor.getAbsoluteColumnIndices(target);
127 }
128 
129 void 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. */
178  _minor.selectFirstRows(k, _container);
179  _minor.selectFirstColumns(k, _container);
180  return true;
181  }
182  else if (_minor.selectNextColumns(k, _container))
183  {
184  /* Here we were able to pick a next subset of columns
185  within the same subset of rows. */
186  return true;
187  }
188  else if (_minor.selectNextRows(k, _container))
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: */
193  _minor.selectFirstColumns(k, _container);
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 
205 bool 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 
218 int 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  int result = 1;
224  for (int k = i - j + 1; k <= i; k++) result *= k;
225  /* Now, result = (i - j + 1) * ... * i. */
226  for (int k = 2; k <= j; k++) result /= k;
227  /* Now, result = (i - j + 1) * ... * i / 1 / 2 ...
228  ... / j = i! / j! / (i - j)!. */
229  return result;
230 }
231 
233 {
234  /* This is a non-recursive implementation. */
235  assume(i >= 0);
236  int result = 1;
237  for (int j = 1; j <= i; j++) result *= j;
238  // Now, result = 1 * 2 * ... * i = i!
239  return result;
240 }
241 
242 int MinorProcessor::NumberOfRetrievals (const int rows, const int columns,
243  const int containerMinorSize,
244  const int minorSize,
245  const bool multipleMinors)
246 {
247  /* This method computes the number of potential retrievals
248  of a single minor when computing all minors of a given size
249  within a matrix of given size. */
250  int result = 0;
251  if (multipleMinors)
252  {
253  /* Here, we would like to compute all minors of size
254  containerMinorSize x containerMinorSize in a matrix
255  of size rows x columns.
256  Then, we need to retrieve any minor of size
257  minorSize x minorSize exactly n times, where n is as
258  follows: */
259  result = IOverJ(rows - minorSize, containerMinorSize - minorSize)
260  * IOverJ(columns - minorSize, containerMinorSize - minorSize)
261  * Faculty(containerMinorSize - minorSize);
262  }
263  else
264  {
265  /* Here, we would like to compute just one minor of size
266  containerMinorSize x containerMinorSize. Then, we need
267  to retrieve any minor of size minorSize x minorSize exactly
268  (containerMinorSize - minorSize)! times: */
269  result = Faculty(containerMinorSize - minorSize);
270  }
271  return result;
272 }
273 
275  _container(0, NULL, 0, NULL),
276  _minor(0, NULL, 0, NULL),
277  _containerRows(0),
278  _containerColumns(0),
279  _minorSize(0),
280  _rows(0),
281  _columns(0)
282 {
283 }
284 
286 
288 {
289  _intMatrix = 0;
290 }
291 
293 {
294  char h[32];
295  string t = "";
296  string s = "IntMinorProcessor:";
297  s += "\n matrix: ";
298  sprintf(h, "%d", _rows); s += h;
299  s += " x ";
300  sprintf(h, "%d", _columns); s += h;
301  for (int r = 0; r < _rows; r++)
302  {
303  s += "\n ";
304  for (int c = 0; c < _columns; c++)
305  {
306  sprintf(h, "%d", getEntry(r, c)); t = h;
307  for (int k = 0; k < int(4 - strlen(h)); k++) s += " ";
308  s += t;
309  }
310  }
311  int myIndexArray[500];
312  s += "\n considered submatrix has row indices: ";
313  _container.getAbsoluteRowIndices(myIndexArray);
314  for (int k = 0; k < _containerRows; k++)
315  {
316  if (k != 0) s += ", ";
317  sprintf(h, "%d", myIndexArray[k]); s += h;
318  }
319  s += " (first row of matrix has index 0)";
320  s += "\n considered submatrix has column indices: ";
321  _container.getAbsoluteColumnIndices(myIndexArray);
322  for (int k = 0; k < _containerColumns; k++)
323  {
324  if (k != 0) s += ", ";
325  sprintf(h, "%d", myIndexArray[k]); s += h;
326  }
327  s += " (first column of matrix has index 0)";
328  s += "\n size of considered minor(s): ";
329  sprintf(h, "%d", _minorSize); s += h;
330  s += "x";
331  s += h;
332  return s;
333 }
334 
336 {
337  /* free memory of _intMatrix */
338  delete [] _intMatrix; _intMatrix = 0;
339 }
340 
341 bool IntMinorProcessor::isEntryZero (const int absoluteRowIndex,
342  const int absoluteColumnIndex) const
343 {
344  return getEntry(absoluteRowIndex, absoluteColumnIndex) == 0;
345 }
346 
347 void IntMinorProcessor::defineMatrix (const int numberOfRows,
348  const int numberOfColumns,
349  const int* matrix)
350 {
351  /* free memory of _intMatrix */
353 
354  _rows = numberOfRows;
355  _columns = numberOfColumns;
356 
357  /* allocate memory for new entries in _intMatrix */
358  int n = _rows * _columns;
359  _intMatrix = (int*)omAlloc(n*sizeof(int));
360 
361  /* copying values from one-dimensional method
362  parameter "matrix" */
363  for (int i = 0; i < n; i++)
364  _intMatrix[i] = matrix[i];
365 }
366 
368  const int* rowIndices,
369  const int* columnIndices,
371  const int characteristic,
372  const ideal& iSB)
373 {
374  defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
376  /* call a helper method which recursively computes the minor using the
377  cache c: */
378  return getMinorPrivateLaplace(dimension, _container, false, c,
379  characteristic, iSB);
380 }
381 
383  const int* rowIndices,
384  const int* columnIndices,
385  const int characteristic,
386  const ideal& iSB,
387  const char* algorithm)
388 {
389  defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
391 
392  /* call a helper method which computes the minor (without a cache): */
393  if (strcmp(algorithm, "Laplace") == 0)
394  return getMinorPrivateLaplace(_minorSize, _container, characteristic,
395  iSB);
396  else if (strcmp(algorithm, "Bareiss") == 0)
397  return getMinorPrivateBareiss(_minorSize, _container, characteristic,
398  iSB);
399  else assume(false);
400 
401  /* The following code is never reached and just there to make the
402  compiler happy: */
403  return IntMinorValue();
404 }
405 
407  const ideal& iSB,
408  const char* algorithm)
409 {
410  /* call a helper method which computes the minor (without a cache): */
411  if (strcmp(algorithm, "Laplace") == 0)
412  return getMinorPrivateLaplace(_minorSize, _minor, characteristic, iSB);
413  else if (strcmp(algorithm, "Bareiss") == 0)
414  return getMinorPrivateBareiss(_minorSize, _minor, characteristic, iSB);
415  else assume(false);
416 
417  /* The following code is never reached and just there to make the
418  compiler happy: */
419  return IntMinorValue();
420 }
421 
423  IntMinorValue>& c,
424  const int characteristic,
425  const ideal& iSB)
426 {
427  /* computation with cache */
428  return getMinorPrivateLaplace(_minorSize, _minor, true, c, characteristic,
429  iSB);
430 }
431 
432 /* computes the reduction of an integer i modulo an ideal
433  which captures a std basis */
434 int getReduction (const int i, const ideal& iSB)
435 {
436  if (i == 0) return 0;
437  poly f = pISet(i);
438  poly g = kNF(iSB, currRing->qideal, f);
439  int result = 0;
440  if (g != NULL) result = n_Int(pGetCoeff(g), currRing->cf);
441  pDelete(&f);
442  pDelete(&g);
443  return result;
444 }
445 
447  const int k,
448  const MinorKey& mk,
449  const int characteristic,
450  const ideal& iSB)
451 {
452  assume(k > 0); /* k is the minor's dimension; the minor must be at least
453  1x1 */
454  /* The method works by recursion, and using Lapace's Theorem along the
455  row/column with the most zeros. */
456  if (k == 1)
457  {
459  if (characteristic != 0) e = e % characteristic;
460  if (iSB != 0) e = getReduction(e, iSB);
461  return IntMinorValue(e, 0, 0, 0, 0, -1, -1); /* "-1" is to signal that any
462  statistics about the number
463  of retrievals does not make
464  sense, as we do not use a
465  cache. */
466  }
467  else
468  {
469  /* Here, the minor must be 2x2 or larger. */
470  int b = getBestLine(k, mk); /* row or column with most
471  zeros */
472  int result = 0; /* This will contain the
473  value of the minor. */
474  int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions and
475  multiplications, ..."a*"
476  for accumulated operation
477  counters */
478  bool hadNonZeroEntry = false;
479  if (b >= 0)
480  {
481  /* This means that the best line is the row with absolute (0-based)
482  index b.
483  Using Laplace, the sign of the contributing minors must be iterating;
484  the initial sign depends on the relative index of b in minorRowKey: */
485  int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
486  for (int c = 0; c < k; c++) /* This iterates over all involved columns. */
487  {
488  int absoluteC = mk.getAbsoluteColumnIndex(c);
489  if (getEntry(b, absoluteC) != 0) /* Only then do we have to consider
490  this sub-determinante. */
491  {
492  hadNonZeroEntry = true;
493  /* Next MinorKey is mk with row b and column absoluteC omitted: */
494  MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
495  IntMinorValue mv = getMinorPrivateLaplace(k - 1, subMk,
496  characteristic, iSB); /* recursive call */
497  m += mv.getMultiplications();
498  s += mv.getAdditions();
500  as += mv.getAccumulatedAdditions();
501  /* adding sub-determinante times matrix entry
502  times appropriate sign: */
503  result += sign * mv.getResult() * getEntry(b, absoluteC);
504 
505  if (characteristic != 0) result = result % characteristic;
506  s++; m++; as++, am++; /* This is for the last addition and
507  multiplication. */
508  }
509  sign = - sign; /* alternating the sign */
510  }
511  }
512  else
513  {
514  b = - b - 1;
515  /* This means that the best line is the column with absolute (0-based)
516  index b.
517  Using Laplace, the sign of the contributing minors must be iterating;
518  the initial sign depends on the relative index of b in
519  minorColumnKey: */
520  int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
521  for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
522  {
523  int absoluteR = mk.getAbsoluteRowIndex(r);
524  if (getEntry(absoluteR, b) != 0) /* Only then do we have to consider
525  this sub-determinante. */
526  {
527  hadNonZeroEntry = true;
528  /* Next MinorKey is mk with row absoluteR and column b omitted. */
529  MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
530  IntMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, characteristic, iSB); /* recursive call */
531  m += mv.getMultiplications();
532  s += mv.getAdditions();
534  as += mv.getAccumulatedAdditions();
535  /* adding sub-determinante times matrix entry
536  times appropriate sign: */
537  result += sign * mv.getResult() * getEntry(absoluteR, b);
538  if (characteristic != 0) result = result % characteristic;
539  s++; m++; as++, am++; /* This is for the last addition and
540  multiplication. */
541  }
542  sign = - sign; /* alternating the sign */
543  }
544  }
545  if (hadNonZeroEntry)
546  {
547  s--; as--; /* first addition was 0 + ..., so we do not count it */
548  }
549  if (s < 0) s = 0; /* may happen when all subminors are zero and no
550  addition needs to be performed */
551  if (as < 0) as = 0; /* may happen when all subminors are zero and no
552  addition needs to be performed */
553  if (iSB != 0) result = getReduction(result, iSB);
554  IntMinorValue newMV(result, m, s, am, as, -1, -1);
555  /* "-1" is to signal that any statistics about the number of retrievals
556  does not make sense, as we do not use a cache. */
557  return newMV;
558  }
559 }
560 
561 /* This method can only be used in the case of coefficients
562  coming from a field or at least from an integral domain. */
564  const int k,
565  const MinorKey& mk,
566  const int characteristic,
567  const ideal& iSB)
568 {
569  assume(k > 0); /* k is the minor's dimension; the minor must be at least
570  1x1 */
571  int *theRows=(int*)omAlloc(k*sizeof(int));
572  mk.getAbsoluteRowIndices(theRows);
573  int *theColumns=(int*)omAlloc(k*sizeof(int));
574  mk.getAbsoluteColumnIndices(theColumns);
575  /* the next line provides the return value for the case k = 1 */
576  int e = getEntry(theRows[0], theColumns[0]);
577  if (characteristic != 0) e = e % characteristic;
578  if (iSB != 0) e = getReduction(e, iSB);
579  IntMinorValue mv(e, 0, 0, 0, 0, -1, -1);
580  if (k > 1)
581  {
582  /* the matrix to perform Bareiss with */
583  long *tempMatrix=(long*)omAlloc(k * k *sizeof(long));
584  /* copy correct set of entries from _intMatrix to tempMatrix */
585  int i = 0;
586  for (int r = 0; r < k; r++)
587  for (int c = 0; c < k; c++)
588  {
589  e = getEntry(theRows[r], theColumns[c]);
590  if (characteristic != 0) e = e % characteristic;
591  tempMatrix[i++] = e;
592  }
593  /* Bareiss algorithm operating on tempMatrix which is at least 2x2 */
594  int sign = 1; /* This will store the correct sign resulting
595  from permuting the rows of tempMatrix. */
596  int *rowPermutation=(int*)omAlloc(k*sizeof(int));
597  /* This is for storing the permutation of rows
598  resulting from searching for a non-zero
599  pivot element. */
600  for (int i = 0; i < k; i++) rowPermutation[i] = i;
601  int divisor = 1; /* the Bareiss divisor */
602  for (int r = 0; r <= k - 2; r++)
603  {
604  /* look for a non-zero entry in column r: */
605  int i = r;
606  while ((i < k) && (tempMatrix[rowPermutation[i] * k + r] == 0))
607  i++;
608  if (i == k)
609  /* There is no non-zero entry; hence the minor is zero. */
610  return IntMinorValue(0, 0, 0, 0, 0, -1, -1);
611  if (i != r)
612  {
613  /* We swap the rows with indices r and i: */
614  int j = rowPermutation[i];
615  rowPermutation[i] = rowPermutation[r];
616  rowPermutation[r] = j;
617  /* Now we know that tempMatrix[rowPermutation[r] * k + r] is not zero.
618  But carefull; we have to negate the sign, as there is always an odd
619  number of row transpositions to swap two given rows of a matrix. */
620  sign = -sign;
621  }
622  if (r >= 1) divisor = tempMatrix[rowPermutation[r - 1] * k + r - 1];
623  for (int rr = r + 1; rr < k; rr++)
624  for (int cc = r + 1; cc < k; cc++)
625  {
626  e = rowPermutation[rr] * k + cc;
627  /* Attention: The following may cause an overflow and
628  thus a wrong result: */
629  tempMatrix[e] = tempMatrix[e] * tempMatrix[rowPermutation[r] * k + r]
630  - tempMatrix[rowPermutation[r] * k + cc]
631  * tempMatrix[rowPermutation[rr] * k + r];
632  /* The following is, by theory, always a division without
633  remainder: */
634  tempMatrix[e] = tempMatrix[e] / divisor;
635  if (characteristic != 0)
636  tempMatrix[e] = tempMatrix[e] % characteristic;
637  }
638  omFree(rowPermutation);
639  omFree(tempMatrix);
640  }
641  int theValue = sign * tempMatrix[rowPermutation[k - 1] * k + k - 1];
642  if (iSB != 0) theValue = getReduction(theValue, iSB);
643  mv = IntMinorValue(theValue, 0, 0, 0, 0, -1, -1);
644  }
645  omFree(theRows);
646  omFree(theColumns);
647  return mv;
648 }
649 
650 int IntMinorProcessor::getEntry (const int rowIndex,
651  const int columnIndex) const
652 {
653  return _intMatrix[rowIndex * _columns + columnIndex];
654 }
655 
657  const int k, const MinorKey& mk,
658  const bool multipleMinors,
660  const int characteristic, const ideal& iSB)
661 {
662  assume(k > 0); /* k is the minor's dimension; the minor must be at least
663  1x1 */
664  /* The method works by recursion, and using Lapace's Theorem along
665  the row/column with the most zeros. */
666  if (k == 1)
667  {
669  if (characteristic != 0) e = e % characteristic;
670  if (iSB != 0) e = getReduction(e, iSB);
671  return IntMinorValue(e, 0, 0, 0, 0, -1, -1);
672  /* we set "-1" as, for k == 1, we do not have any cache retrievals */
673  }
674  else
675  {
676  int b = getBestLine(k, mk); /* row or column with
677  most zeros */
678  int result = 0; /* This will contain the
679  value of the minor. */
680  int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions
681  and multiplications,
682  ..."a*" for
683  accumulated operation
684  counters */
685  IntMinorValue mv(0, 0, 0, 0, 0, 0, 0); /* for storing all
686  intermediate minors */
687  bool hadNonZeroEntry = false;
688  if (b >= 0)
689  {
690  /* This means that the best line is the row with absolute (0-based)
691  index b.
692  Using Laplace, the sign of the contributing minors must be
693  iterating; the initial sign depends on the relative index of b
694  in minorRowKey: */
695  int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
696  for (int c = 0; c < k; c++) /* This iterates over all involved
697  columns. */
698  {
699  int absoluteC = mk.getAbsoluteColumnIndex(c);
700  if (getEntry(b, absoluteC) != 0) /* Only then do we have to consider
701  this sub-determinante. */
702  {
703  hadNonZeroEntry = true;
704  /* Next MinorKey is mk with row b and column absoluteC omitted. */
705  MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
706  if (cch.hasKey(subMk))
707  { /* trying to find the result in the cache */
708  mv = cch.getValue(subMk);
709  mv.incrementRetrievals(); /* once more, we made use of the cached
710  value for key mk */
711  cch.put(subMk, mv); /* We need to do this with "put", as the
712  (altered) number of retrievals may have
713  an impact on the internal ordering among
714  the cached entries. */
715  }
716  else
717  {
718  mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch,
719  characteristic, iSB); /* recursive call */
720  /* As this minor was not in the cache, we count the additions
721  and multiplications that we needed to perform in the
722  recursive call: */
723  m += mv.getMultiplications();
724  s += mv.getAdditions();
725  }
726  /* In any case, we count all nested operations in the accumulative
727  counters: */
729  as += mv.getAccumulatedAdditions();
730  /* adding sub-determinante times matrix entry times appropriate
731  sign */
732  result += sign * mv.getResult() * getEntry(b, absoluteC);
733  if (characteristic != 0) result = result % characteristic;
734  s++; m++; as++; am++; /* This is for the last addition and
735  multiplication. */
736  }
737  sign = - sign; /* alternating the sign */
738  }
739  }
740  else
741  {
742  b = - b - 1;
743  /* This means that the best line is the column with absolute (0-based)
744  index b.
745  Using Laplace, the sign of the contributing minors must be iterating;
746  the initial sign depends on the relative index of b in
747  minorColumnKey: */
748  int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
749  for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
750  {
751  int absoluteR = mk.getAbsoluteRowIndex(r);
752  if (getEntry(absoluteR, b) != 0) /* Only then do we have to consider
753  this sub-determinante. */
754  {
755  hadNonZeroEntry = true;
756  /* Next MinorKey is mk with row absoluteR and column b omitted. */
757  MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
758  if (cch.hasKey(subMk))
759  { /* trying to find the result in the cache */
760  mv = cch.getValue(subMk);
761  mv.incrementRetrievals(); /* once more, we made use of the cached
762  value for key mk */
763  cch.put(subMk, mv); /* We need to do this with "put", as the
764  (altered) number of retrievals may have an
765  impact on the internal ordering among the
766  cached entries. */
767  }
768  else
769  {
770  mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch, characteristic, iSB); /* recursive call */
771  /* As this minor was not in the cache, we count the additions and
772  multiplications that we needed to do in the recursive call: */
773  m += mv.getMultiplications();
774  s += mv.getAdditions();
775  }
776  /* In any case, we count all nested operations in the accumulative
777  counters: */
779  as += mv.getAccumulatedAdditions();
780  /* adding sub-determinante times matrix entry times appropriate
781  sign: */
782  result += sign * mv.getResult() * getEntry(absoluteR, b);
783  if (characteristic != 0) result = result % characteristic;
784  s++; m++; as++; am++; /* This is for the last addition and
785  multiplication. */
786  }
787  sign = - sign; /* alternating the sign */
788  }
789  }
790  /* Let's cache the newly computed minor: */
791  int potentialRetrievals = NumberOfRetrievals(_containerRows,
793  _minorSize, k,
794  multipleMinors);
795  if (hadNonZeroEntry)
796  {
797  s--; as--; /* first addition was 0 + ..., so we do not count it */
798  }
799  if (s < 0) s = 0; /* may happen when all subminors are zero and no
800  addition needs to be performed */
801  if (as < 0) as = 0; /* may happen when all subminors are zero and no
802  addition needs to be performed */
803  if (iSB != 0) result = getReduction(result, iSB);
804  IntMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
805  cch.put(mk, newMV); /* Here's the actual put inside the cache. */
806  return newMV;
807  }
808 }
809 
811 {
812  _polyMatrix = NULL;
813 }
814 
815 poly PolyMinorProcessor::getEntry (const int rowIndex,
816  const int columnIndex) const
817 {
818  return _polyMatrix[rowIndex * _columns + columnIndex];
819 }
820 
821 bool PolyMinorProcessor::isEntryZero (const int absoluteRowIndex,
822  const int absoluteColumnIndex) const
823 {
824  return getEntry(absoluteRowIndex, absoluteColumnIndex) == NULL;
825 }
826 
828 {
829  char h[32];
830  string t = "";
831  string s = "PolyMinorProcessor:";
832  s += "\n matrix: ";
833  sprintf(h, "%d", _rows); s += h;
834  s += " x ";
835  sprintf(h, "%d", _columns); s += h;
836  int myIndexArray[500];
837  s += "\n considered submatrix has row indices: ";
838  _container.getAbsoluteRowIndices(myIndexArray);
839  for (int k = 0; k < _containerRows; k++)
840  {
841  if (k != 0) s += ", ";
842  sprintf(h, "%d", myIndexArray[k]); s += h;
843  }
844  s += " (first row of matrix has index 0)";
845  s += "\n considered submatrix has column indices: ";
846  _container.getAbsoluteColumnIndices(myIndexArray);
847  for (int k = 0; k < _containerColumns; k++)
848  {
849  if (k != 0) s += ", ";
850  sprintf(h, "%d", myIndexArray[k]); s += h;
851  }
852  s += " (first column of matrix has index 0)";
853  s += "\n size of considered minor(s): ";
854  sprintf(h, "%d", _minorSize); s += h;
855  s += "x";
856  s += h;
857  return s;
858 }
859 
861 {
862  /* free memory of _polyMatrix */
863  int n = _rows * _columns;
864  for (int i = 0; i < n; i++)
867 }
868 
869 void PolyMinorProcessor::defineMatrix (const int numberOfRows,
870  const int numberOfColumns,
871  const poly* polyMatrix)
872 {
873  /* free memory of _polyMatrix */
874  int n = _rows * _columns;
875  for (int i = 0; i < n; i++)
878 
879  _rows = numberOfRows;
880  _columns = numberOfColumns;
881  n = _rows * _columns;
882 
883  /* allocate memory for new entries in _polyMatrix */
884  _polyMatrix = (poly*)omAlloc(n*sizeof(poly));
885 
886  /* copying values from one-dimensional method
887  parameter "polyMatrix" */
888  for (int i = 0; i < n; i++)
889  _polyMatrix[i] = pCopy(polyMatrix[i]);
890 }
891 
893  const int* rowIndices,
894  const int* columnIndices,
896  const ideal& iSB)
897 {
898  defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
900  /* call a helper method which recursively computes the minor using the cache
901  c: */
902  return getMinorPrivateLaplace(dimension, _container, false, c, iSB);
903 }
904 
906  const int* rowIndices,
907  const int* columnIndices,
908  const char* algorithm,
909  const ideal& iSB)
910 {
911  defineSubMatrix(dimension, rowIndices, dimension, columnIndices);
913  /* call a helper method which computes the minor (without using a cache): */
914  if (strcmp(algorithm, "Laplace") == 0)
916  else if (strcmp(algorithm, "Bareiss") == 0)
918  else assume(false);
919 
920  /* The following code is never reached and just there to make the
921  compiler happy: */
922  return PolyMinorValue();
923 }
924 
926  const ideal& iSB)
927 {
928  /* call a helper method which computes the minor (without using a
929  cache): */
930  if (strcmp(algorithm, "Laplace") == 0)
932  else if (strcmp(algorithm, "Bareiss") == 0)
934  else assume(false);
935 
936  /* The following code is never reached and just there to make the
937  compiler happy: */
938  return PolyMinorValue();
939 }
940 
942  PolyMinorValue>& c,
943  const ideal& iSB)
944 {
945  /* computation with cache */
946  return getMinorPrivateLaplace(_minorSize, _minor, true, c, iSB);
947 }
948 
949 /* assumes that all entries in polyMatrix are already reduced w.r.t. iSB */
951  const MinorKey& mk,
952  const ideal& iSB)
953 {
954  assume(k > 0); /* k is the minor's dimension; the minor must be at least
955  1x1 */
956  /* The method works by recursion, and using Lapace's Theorem along the
957  row/column with the most zeros. */
958  if (k == 1)
959  {
961  mk.getAbsoluteColumnIndex(0)),
962  0, 0, 0, 0, -1, -1);
963  /* "-1" is to signal that any statistics about the number of retrievals
964  does not make sense, as we do not use a cache. */
965  return pmv;
966  }
967  else
968  {
969  /* Here, the minor must be 2x2 or larger. */
970  int b = getBestLine(k, mk); /* row or column with most
971  zeros */
972  poly result = NULL; /* This will contain the
973  value of the minor. */
974  int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions
975  and multiplications,
976  ..."a*" for accumulated
977  operation counters */
978  bool hadNonZeroEntry = false;
979  if (b >= 0)
980  {
981  /* This means that the best line is the row with absolute (0-based)
982  index b.
983  Using Laplace, the sign of the contributing minors must be iterating;
984  the initial sign depends on the relative index of b in minorRowKey: */
985  int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
986  for (int c = 0; c < k; c++) /* This iterates over all involved columns. */
987  {
988  int absoluteC = mk.getAbsoluteColumnIndex(c);
989  if (!isEntryZero(b, absoluteC)) /* Only then do we have to consider
990  this sub-determinante. */
991  {
992  hadNonZeroEntry = true;
993  /* Next MinorKey is mk with row b and column absoluteC omitted. */
994  MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
995  /* recursive call: */
996  PolyMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, iSB);
997  m += mv.getMultiplications();
998  s += mv.getAdditions();
1000  as += mv.getAccumulatedAdditions();
1001  poly signPoly = pISet(sign);
1002  poly temp = pp_Mult_qq(mv.getResult(), getEntry(b, absoluteC),
1003  currRing);
1004  temp = p_Mult_q(signPoly, temp, currRing);
1005  result = p_Add_q(result, temp, currRing);
1006 #ifdef COUNT_AND_PRINT_OPERATIONS
1007  multsPoly++;
1008  addsPoly++;
1009  multsMon += pLength(mv.getResult()) * pLength(getEntry(b, absoluteC));
1010 #endif
1011  s++; m++; as++, am++; /* This is for the addition and multiplication
1012  in the previous lines of code. */
1013  }
1014  sign = - sign; /* alternating the sign */
1015  }
1016  }
1017  else
1018  {
1019  b = - b - 1;
1020  /* This means that the best line is the column with absolute (0-based)
1021  index b.
1022  Using Laplace, the sign of the contributing minors must be iterating;
1023  the initial sign depends on the relative index of b in
1024  minorColumnKey: */
1025  int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
1026  for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
1027  {
1028  int absoluteR = mk.getAbsoluteRowIndex(r);
1029  if (!isEntryZero(absoluteR, b)) /* Only then do we have to consider
1030  this sub-determinante. */
1031  {
1032  hadNonZeroEntry = true;
1033  /* This is mk with row absoluteR and column b omitted. */
1034  MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
1035  /* recursive call: */
1036  PolyMinorValue mv = getMinorPrivateLaplace(k - 1, subMk, iSB);
1037  m += mv.getMultiplications();
1038  s += mv.getAdditions();
1039  am += mv.getAccumulatedMultiplications();
1040  as += mv.getAccumulatedAdditions();
1041  poly signPoly = pISet(sign);
1042  poly temp = pp_Mult_qq(mv.getResult(), getEntry(absoluteR, b),
1043  currRing);
1044  temp = p_Mult_q(signPoly, temp, currRing);
1045  result = p_Add_q(result, temp, currRing);
1046 #ifdef COUNT_AND_PRINT_OPERATIONS
1047  multsPoly++;
1048  addsPoly++;
1049  multsMon += pLength(mv.getResult()) * pLength(getEntry(absoluteR, b));
1050 #endif
1051  s++; m++; as++, am++; /* This is for the addition and multiplication
1052  in the previous lines of code. */
1053  }
1054  sign = - sign; /* alternating the sign */
1055  }
1056  }
1057  if (hadNonZeroEntry)
1058  {
1059  s--; as--; /* first addition was 0 + ..., so we do not count it */
1060 #ifdef COUNT_AND_PRINT_OPERATIONS
1061  addsPoly--;
1062 #endif
1063  }
1064  if (s < 0) s = 0; /* may happen when all subminors are zero and no
1065  addition needs to be performed */
1066  if (as < 0) as = 0; /* may happen when all subminors are zero and no
1067  addition needs to be performed */
1068  if (iSB != NULL)
1069  {
1070  poly tmpresult = kNF(iSB, currRing->qideal, result);
1071  pDelete(&result);
1072  result=tmpresult;
1073  }
1074  PolyMinorValue newMV(result, m, s, am, as, -1, -1);
1075  /* "-1" is to signal that any statistics about the number of retrievals
1076  does not make sense, as we do not use a cache. */
1077  pDelete(&result);
1078  return newMV;
1079  }
1080 }
1081 
1082 /* assumes that all entries in polyMatrix are already reduced w.r.t. iSB */
1084  const int k,
1085  const MinorKey& mk,
1086  const bool multipleMinors,
1088  const ideal& iSB)
1089 {
1090  assume(k > 0); /* k is the minor's dimension; the minor must be at least
1091  1x1 */
1092  /* The method works by recursion, and using Lapace's Theorem along
1093  the row/column with the most zeros. */
1094  if (k == 1)
1095  {
1097  mk.getAbsoluteColumnIndex(0)),
1098  0, 0, 0, 0, -1, -1);
1099  /* we set "-1" as, for k == 1, we do not have any cache retrievals */
1100  return pmv;
1101  }
1102  else
1103  {
1104  int b = getBestLine(k, mk); /* row or column with most
1105  zeros */
1106  poly result = NULL; /* This will contain the
1107  value of the minor. */
1108  int s = 0; int m = 0; int as = 0; int am = 0; /* counters for additions
1109  and multiplications,
1110  ..."a*" for accumulated
1111  operation counters */
1112  bool hadNonZeroEntry = false;
1113  if (b >= 0)
1114  {
1115  /* This means that the best line is the row with absolute (0-based)
1116  index b.
1117  Using Laplace, the sign of the contributing minors must be iterating;
1118  the initial sign depends on the relative index of b in
1119  minorRowKey: */
1120  int sign = (mk.getRelativeRowIndex(b) % 2 == 0 ? 1 : -1);
1121  for (int c = 0; c < k; c++) /* This iterates over all involved columns. */
1122  {
1123  int absoluteC = mk.getAbsoluteColumnIndex(c);
1124  if (!isEntryZero(b, absoluteC)) /* Only then do we have to consider
1125  this sub-determinante. */
1126  {
1127  hadNonZeroEntry = true;
1128  PolyMinorValue mv; /* for storing all intermediate minors */
1129  /* Next MinorKey is mk with row b and column absoluteC omitted. */
1130  MinorKey subMk = mk.getSubMinorKey(b, absoluteC);
1131  if (cch.hasKey(subMk))
1132  { /* trying to find the result in the cache */
1133  mv = cch.getValue(subMk);
1134  mv.incrementRetrievals(); /* once more, we made use of the cached
1135  value for key mk */
1136  cch.put(subMk, mv); /* We need to do this with "put", as the
1137  (altered) number of retrievals may have an
1138  impact on the internal ordering among cache
1139  entries. */
1140  }
1141  else
1142  {
1143  /* recursive call: */
1144  mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch,
1145  iSB);
1146  /* As this minor was not in the cache, we count the additions and
1147  multiplications that we needed to do in the recursive call: */
1148  m += mv.getMultiplications();
1149  s += mv.getAdditions();
1150  }
1151  /* In any case, we count all nested operations in the accumulative
1152  counters: */
1153  am += mv.getAccumulatedMultiplications();
1154  as += mv.getAccumulatedAdditions();
1155  poly signPoly = pISet(sign);
1156  poly temp = pp_Mult_qq(mv.getResult(), getEntry(b, absoluteC),
1157  currRing);
1158  temp = p_Mult_q(signPoly, temp, currRing);
1159  result = p_Add_q(result, temp, currRing);
1160 #ifdef COUNT_AND_PRINT_OPERATIONS
1161  multsPoly++;
1162  addsPoly++;
1163  multsMon += pLength(mv.getResult()) * pLength(getEntry(b, absoluteC));
1164 #endif
1165  s++; m++; as++; am++; /* This is for the addition and multiplication
1166  in the previous lines of code. */
1167  }
1168  sign = - sign; /* alternating the sign */
1169  }
1170  }
1171  else
1172  {
1173  b = - b - 1;
1174  /* This means that the best line is the column with absolute (0-based)
1175  index b.
1176  Using Laplace, the sign of the contributing minors must be iterating;
1177  the initial sign depends on the relative index of b in
1178  minorColumnKey: */
1179  int sign = (mk.getRelativeColumnIndex(b) % 2 == 0 ? 1 : -1);
1180  for (int r = 0; r < k; r++) /* This iterates over all involved rows. */
1181  {
1182  int absoluteR = mk.getAbsoluteRowIndex(r);
1183  if (!isEntryZero(absoluteR, b)) /* Only then do we have to consider
1184  this sub-determinante. */
1185  {
1186  hadNonZeroEntry = true;
1187  PolyMinorValue mv; /* for storing all intermediate minors */
1188  /* Next MinorKey is mk with row absoluteR and column b omitted. */
1189  MinorKey subMk = mk.getSubMinorKey(absoluteR, b);
1190  if (cch.hasKey(subMk))
1191  { /* trying to find the result in the cache */
1192  mv = cch.getValue(subMk);
1193  mv.incrementRetrievals(); /* once more, we made use of the cached
1194  value for key mk */
1195  cch.put(subMk, mv); /* We need to do this with "put", as the
1196  (altered) number of retrievals may have an
1197  impact on the internal ordering among the
1198  cached entries. */
1199  }
1200  else
1201  {
1202  mv = getMinorPrivateLaplace(k - 1, subMk, multipleMinors, cch,
1203  iSB); /* recursive call */
1204  /* As this minor was not in the cache, we count the additions and
1205  multiplications that we needed to do in the recursive call: */
1206  m += mv.getMultiplications();
1207  s += mv.getAdditions();
1208  }
1209  /* In any case, we count all nested operations in the accumulative
1210  counters: */
1211  am += mv.getAccumulatedMultiplications();
1212  as += mv.getAccumulatedAdditions();
1213  poly signPoly = pISet(sign);
1214  poly temp = pp_Mult_qq(mv.getResult(), getEntry(absoluteR, b),
1215  currRing);
1216  temp = p_Mult_q(signPoly, temp, currRing);
1217  result = p_Add_q(result, temp, currRing);
1218 #ifdef COUNT_AND_PRINT_OPERATIONS
1219  multsPoly++;
1220  addsPoly++;
1221  multsMon += pLength(mv.getResult()) * pLength(getEntry(absoluteR, b));
1222 #endif
1223  s++; m++; as++; am++; /* This is for the addition and multiplication
1224  in the previous lines of code. */
1225  }
1226  sign = - sign; /* alternating the sign */
1227  }
1228  }
1229  /* Let's cache the newly computed minor: */
1230  int potentialRetrievals = NumberOfRetrievals(_containerRows,
1232  _minorSize,
1233  k,
1234  multipleMinors);
1235  if (hadNonZeroEntry)
1236  {
1237  s--; as--; /* first addition was 0 + ..., so we do not count it */
1238 #ifdef COUNT_AND_PRINT_OPERATIONS
1239  addsPoly--;
1240 #endif
1241  }
1242  if (s < 0) s = 0; /* may happen when all subminors are zero and no
1243  addition needs to be performed */
1244  if (as < 0) as = 0; /* may happen when all subminors are zero and no
1245  addition needs to be performed */
1246  if (iSB != NULL)
1247  {
1248  poly tmpresult = kNF(iSB, currRing->qideal, result);
1249  pDelete(&tmpresult);
1250  result=tmpresult;
1251  }
1252  PolyMinorValue newMV(result, m, s, am, as, 1, potentialRetrievals);
1253  pDelete(&result); result = NULL;
1254  cch.put(mk, newMV); /* Here's the actual put inside the cache. */
1255  return newMV;
1256  }
1257 }
1258 
1259 /* This can only be used in the case of coefficients coming from a field
1260  or at least an integral domain. */
1261 static void addOperationBucket(poly f1, poly f2, kBucket_pt bucket)
1262 {
1263  /* fills all terms of f1 * f2 into the bucket */
1264  poly a = f1; poly b = f2;
1265  int aLen = pLength(a); int bLen = pLength(b);
1266  if (aLen > bLen)
1267  {
1268  b = f1; a = f2; bLen = aLen;
1269  }
1270  pNormalize(b);
1271 
1272  while (a != NULL)
1273  {
1274  /* The next line actually uses only LT(a): */
1275  kBucket_Plus_mm_Mult_pp(bucket, a, b, bLen);
1276  a = pNext(a);
1277  }
1278 }
1279 
1280 /* computes the polynomial (p1 * p2 - p3 * p4) and puts result into p1;
1281  the method destroys the old value of p1;
1282  p2, p3, and p4 may be pNormalize-d but must, apart from that,
1283  not be changed;
1284  This can only be used in the case of coefficients coming from a field
1285  or at least an integral domain. */
1286 static void elimOperationBucketNoDiv(poly &p1, poly p2, poly p3, poly p4)
1287 {
1288 #ifdef COUNT_AND_PRINT_OPERATIONS
1289  if ((pLength(p1) != 0) && (pLength(p2) != 0))
1290  {
1291  multsPoly++;
1292  multsMon += pLength(p1) * pLength(p2);
1293  }
1294  if ((pLength(p3) != 0) && (pLength(p4) != 0))
1295  {
1296  multsPoly++;
1297  multsMon += pLength(p3) * pLength(p4);
1298  }
1299  if ((pLength(p1) != 0) && (pLength(p2) != 0) &&
1300  (pLength(p3) != 0) && (pLength(p4) != 0))
1301  addsPoly++;
1302 #endif
1303  kBucket_pt myBucket = kBucketCreate(currRing);
1304  addOperationBucket(p1, p2, myBucket);
1305  poly p3Neg = pNeg(pCopy(p3));
1306  addOperationBucket(p3Neg, p4, myBucket);
1307  pDelete(&p3Neg);
1308  pDelete(&p1);
1309  p1 = kBucketClear(myBucket);
1310  kBucketDestroy(&myBucket);
1311 }
1312 
1313 /* computes the polynomial (p1 * p2 - p3 * p4) / p5 and puts result into p1;
1314  the method destroys the old value of p1;
1315  p2, p3, p4, and p5 may be pNormalize-d but must, apart from that,
1316  not be changed;
1317  c5 is assumed to be the leading coefficient of p5;
1318  p5Len is assumed to be the length of p5;
1319  This can only be used in the case of coefficients coming from a field
1320  or at least an integral domain. */
1321 void elimOperationBucket(poly &p1, poly &p2, poly &p3, poly &p4, poly &p5,
1322  number &c5, int p5Len)
1323 {
1324 #ifdef COUNT_AND_PRINT_OPERATIONS
1325  if ((pLength(p1) != 0) && (pLength(p2) != 0))
1326  {
1327  multsPoly++;
1328  multsMon += pLength(p1) * pLength(p2);
1329  }
1330  if ((pLength(p3) != 0) && (pLength(p4) != 0))
1331  {
1332  multsPoly++;
1333  multsMon += pLength(p3) * pLength(p4);
1334  }
1335  if ((pLength(p1) != 0) && (pLength(p2) != 0) &&
1336  (pLength(p3) != 0) && (pLength(p4) != 0))
1337  addsPoly++;
1338 #endif
1339  kBucket_pt myBucket = kBucketCreate(currRing);
1340  addOperationBucket(p1, p2, myBucket);
1341  poly p3Neg = pNeg(pCopy(p3));
1342  addOperationBucket(p3Neg, p4, myBucket);
1343  pDelete(&p3Neg);
1344 
1345  /* Now, myBucket contains all terms of p1 * p2 - p3 * p4.
1346  Now we need to perform the polynomial division myBucket / p5
1347  which is known to work without remainder: */
1348  pDelete(&p1); poly helperPoly = NULL;
1349 
1350  poly bucketLm = pCopy(kBucketGetLm(myBucket));
1351  while (bucketLm != NULL)
1352  {
1353  /* divide bucketLm by the leading term of p5 and put result into bucketLm;
1354  we start with the coefficients;
1355  note that bucketLm will always represent a term */
1356  number coeff = nDiv(pGetCoeff(bucketLm), c5);
1357  nNormalize(coeff);
1358  pSetCoeff(bucketLm, coeff);
1359  /* subtract exponent vector of p5 from that of quotient; modifies
1360  quotient */
1361  p_ExpVectorSub(bucketLm, p5, currRing);
1362 #ifdef COUNT_AND_PRINT_OPERATIONS
1363  divsMon++;
1364  multsMonForDiv += p5Len;
1365  multsMon += p5Len;
1366  savedMultsMFD++;
1367  multsPoly++;
1368  multsPolyForDiv++;
1369  addsPoly++;
1370  addsPolyForDiv++;
1371 #endif
1372  kBucket_Minus_m_Mult_p(myBucket, bucketLm, p5, &p5Len);
1373  /* The following lines make bucketLm the new leading term of p1,
1374  i.e., put bucketLm in front of everything which is already in p1.
1375  Thus, after the while loop, we need to revert p1. */
1376  helperPoly = bucketLm;
1377  helperPoly->next = p1;
1378  p1 = helperPoly;
1379 
1380  bucketLm = pCopy(kBucketGetLm(myBucket));
1381  }
1382  p1 = pReverse(p1);
1383  kBucketDestroy(&myBucket);
1384 }
1385 
1386 /* assumes that all entries in polyMatrix are already reduced w.r.t. iSB
1387  This can only be used in the case of coefficients coming from a field!!! */
1389  const MinorKey& mk,
1390  const ideal& iSB)
1391 {
1392  assume(k > 0); /* k is the minor's dimension; the minor must be at least
1393  1x1 */
1394  int *theRows=(int*)omAlloc(k*sizeof(int));
1395  mk.getAbsoluteRowIndices(theRows);
1396  int *theColumns=(int*)omAlloc(k*sizeof(int));
1397  mk.getAbsoluteColumnIndices(theColumns);
1398  if (k == 1)
1399  {
1400  PolyMinorValue tmp=PolyMinorValue(getEntry(theRows[0], theColumns[0]),
1401  0, 0, 0, 0, -1, -1);
1402  omFree(theColumns);
1403  omFree(theRows);
1404  return tmp;
1405  }
1406  else /* k > 0 */
1407  {
1408  /* the matrix to perform Bareiss with */
1409  poly* tempMatrix = (poly*)omAlloc(k * k * sizeof(poly));
1410  /* copy correct set of entries from _polyMatrix to tempMatrix */
1411  int i = 0;
1412  for (int r = 0; r < k; r++)
1413  for (int c = 0; c < k; c++)
1414  tempMatrix[i++] = pCopy(getEntry(theRows[r], theColumns[c]));
1415 
1416  /* Bareiss algorithm operating on tempMatrix which is at least 2x2 */
1417  int sign = 1; /* This will store the correct sign resulting from
1418  permuting the rows of tempMatrix. */
1419  int *rowPermutation=(int*)omAlloc(k*sizeof(int));
1420  /* This is for storing the permutation of rows
1421  resulting from searching for a non-zero pivot
1422  element. */
1423  for (int i = 0; i < k; i++) rowPermutation[i] = i;
1424  poly divisor = NULL;
1425  int divisorLength = 0;
1426  number divisorLC;
1427  for (int r = 0; r <= k - 2; r++)
1428  {
1429  /* look for a non-zero entry in column r, rows = r .. (k - 1)
1430  s.t. the polynomial has least complexity: */
1431  int minComplexity = -1; int complexity = 0; int bestRow = -1;
1432  poly pp = NULL;
1433  for (int i = r; i < k; i++)
1434  {
1435  pp = tempMatrix[rowPermutation[i] * k + r];
1436  if (pp != NULL)
1437  {
1438  if (minComplexity == -1)
1439  {
1440  minComplexity = pSize(pp);
1441  bestRow = i;
1442  }
1443  else
1444  {
1445  complexity = 0;
1446  while ((pp != NULL) && (complexity < minComplexity))
1447  {
1448  complexity += nSize(pGetCoeff(pp)); pp = pNext(pp);
1449  }
1450  if (complexity < minComplexity)
1451  {
1452  minComplexity = complexity;
1453  bestRow = i;
1454  }
1455  }
1456  if (minComplexity <= 1) break; /* terminate the search */
1457  }
1458  }
1459  if (bestRow == -1)
1460  {
1461  /* There is no non-zero entry; hence the minor is zero. */
1462  for (int i = 0; i < k * k; i++) pDelete(&tempMatrix[i]);
1463  return PolyMinorValue(NULL, 0, 0, 0, 0, -1, -1);
1464  }
1465  pNormalize(tempMatrix[rowPermutation[bestRow] * k + r]);
1466  if (bestRow != r)
1467  {
1468  /* We swap the rows with indices r and i: */
1469  int j = rowPermutation[bestRow];
1470  rowPermutation[bestRow] = rowPermutation[r];
1471  rowPermutation[r] = j;
1472  /* Now we know that tempMatrix[rowPermutation[r] * k + r] is not zero.
1473  But carefull; we have to negate the sign, as there is always an odd
1474  number of row transpositions to swap two given rows of a matrix. */
1475  sign = -sign;
1476  }
1477 #if (defined COUNT_AND_PRINT_OPERATIONS) && (COUNT_AND_PRINT_OPERATIONS > 2)
1478  poly w = NULL; int wl = 0;
1479  printf("matrix after %d steps:\n", r);
1480  for (int u = 0; u < k; u++)
1481  {
1482  for (int v = 0; v < k; v++)
1483  {
1484  if ((v < r) && (u > v))
1485  wl = 0;
1486  else
1487  {
1488  w = tempMatrix[rowPermutation[u] * k + v];
1489  wl = pLength(w);
1490  }
1491  printf("%5d ", wl);
1492  }
1493  printf("\n");
1494  }
1495  printCounters ("", false);
1496 #endif
1497  if (r != 0)
1498  {
1499  divisor = tempMatrix[rowPermutation[r - 1] * k + r - 1];
1500  pNormalize(divisor);
1501  divisorLength = pLength(divisor);
1502  divisorLC = pGetCoeff(divisor);
1503  }
1504  for (int rr = r + 1; rr < k; rr++)
1505  for (int cc = r + 1; cc < k; cc++)
1506  {
1507  if (r == 0)
1508  elimOperationBucketNoDiv(tempMatrix[rowPermutation[rr] * k + cc],
1509  tempMatrix[rowPermutation[r] * k + r],
1510  tempMatrix[rowPermutation[r] * k + cc],
1511  tempMatrix[rowPermutation[rr] * k + r]);
1512  else
1513  elimOperationBucket(tempMatrix[rowPermutation[rr] * k + cc],
1514  tempMatrix[rowPermutation[r] * k + r],
1515  tempMatrix[rowPermutation[r] * k + cc],
1516  tempMatrix[rowPermutation[rr] * k + r],
1517  divisor, divisorLC, divisorLength);
1518  }
1519  }
1520 #if (defined COUNT_AND_PRINT_OPERATIONS) && (COUNT_AND_PRINT_OPERATIONS > 2)
1521  poly w = NULL; int wl = 0;
1522  printf("matrix after %d steps:\n", k - 1);
1523  for (int u = 0; u < k; u++)
1524  {
1525  for (int v = 0; v < k; v++)
1526  {
1527  if ((v < k - 1) && (u > v))
1528  wl = 0;
1529  else
1530  {
1531  w = tempMatrix[rowPermutation[u] * k + v];
1532  wl = pLength(w);
1533  }
1534  printf("%5d ", wl);
1535  }
1536  printf("\n");
1537  }
1538 #endif
1539  poly result = tempMatrix[rowPermutation[k - 1] * k + k - 1];
1540  tempMatrix[rowPermutation[k - 1] * k + k - 1]=NULL; /*value now in result*/
1541  if (sign == -1) result = pNeg(result);
1542  if (iSB != NULL)
1543  {
1544  poly tmpresult = kNF(iSB, currRing->qideal, result);
1545  pDelete(&result);
1546  result=tmpresult;
1547  }
1548  PolyMinorValue mv(result, 0, 0, 0, 0, -1, -1);
1549  for (int i = 0; i < k * k; i++) pDelete(&tempMatrix[i]);
1550  omFreeSize(tempMatrix, k * k * sizeof(poly));
1551  omFree(rowPermutation);
1552  omFree(theColumns);
1553  omFree(theRows);
1554  return mv;
1555  }
1556 }
1557 
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
std::string toString(const gfan::ZCone *const c)
Definition: bbcone.cc:27
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
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 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
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 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:547
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
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:827
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:3167
#define assume(x)
Definition: mod2.h:387
#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 poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:908
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1086
static void p_ExpVectorSub(poly p1, poly p2, const ring r)
Definition: p_polys.h:1412
static poly pReverse(poly p)
Definition: p_polys.h:335
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:873
static unsigned pLength(poly a)
Definition: p_polys.h:191
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1123
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatiblity 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:3372