source: git/kernel/linear_algebra/MinorProcessor.cc @ a3f0fea

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