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