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