# source:git/Singular/MinorProcessor.cc@b1dfaf

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