source: git/kernel/linear_algebra/Minor.cc @ c858487

spielwiese
Last change on this file since c858487 was bcd835, checked in by Hans Schoenemann <hannes@…>, 3 years ago
fix: minor: use of omFree/omfree
  • Property mode set to 100644
File size: 37.1 KB
Line 
1
2
3
4#include "kernel/mod2.h"
5
6#include "kernel/linear_algebra/Minor.h"
7
8#include "kernel/structs.h"
9#include "kernel/polys.h"
10
11using namespace std;
12
13void MinorKey::reset()
14{
15  _numberOfRowBlocks = 0;
16  _numberOfColumnBlocks = 0;
17  omfree(_rowKey);
18  _rowKey = NULL;
19  omfree(_columnKey);
20  _columnKey = NULL;
21}
22
23MinorKey::MinorKey (const MinorKey& mk)
24{
25  _numberOfRowBlocks = mk.getNumberOfRowBlocks();
26  _numberOfColumnBlocks = mk.getNumberOfColumnBlocks();;
27
28  /* allocate memory for new entries in _rowKey and _columnKey */
29  _rowKey = (unsigned*)omAlloc(_numberOfRowBlocks*sizeof(unsigned));
30  _columnKey = (unsigned*)omAlloc(_numberOfColumnBlocks*sizeof(unsigned));
31
32  /* copying values from parameter arrays to private arrays */
33  for (int r = 0; r < _numberOfRowBlocks; r++)
34    _rowKey[r] = mk.getRowKey(r);
35  for (int c = 0; c < _numberOfColumnBlocks; c++)
36    _columnKey[c] = mk.getColumnKey(c);
37}
38
39MinorKey& MinorKey::operator=(const MinorKey& mk)
40{
41  omfree(_rowKey); _rowKey = NULL;
42  omfree(_columnKey); _columnKey = NULL;
43  _numberOfRowBlocks = 0;
44  _numberOfColumnBlocks = 0;
45
46  _numberOfRowBlocks = mk.getNumberOfRowBlocks();
47  _numberOfColumnBlocks = mk.getNumberOfColumnBlocks();;
48
49  /* allocate memory for new entries in _rowKey and _columnKey */
50  _rowKey = (unsigned*)omalloc(_numberOfRowBlocks*sizeof(unsigned));
51  _columnKey = (unsigned*)omalloc(_numberOfColumnBlocks*sizeof(unsigned));
52
53  /* copying values from parameter arrays to private arrays */
54  for (int r = 0; r < _numberOfRowBlocks; r++)
55      _rowKey[r] = mk.getRowKey(r);
56  for (int c = 0; c < _numberOfColumnBlocks; c++)
57      _columnKey[c] = mk.getColumnKey(c);
58
59  return *this;
60}
61
62void MinorKey::set(const int lengthOfRowArray, const unsigned int* rowKey,
63                   const int lengthOfColumnArray,
64                   const unsigned int* columnKey)
65{
66  /* free memory of _rowKey and _columnKey */
67  if (_numberOfRowBlocks > 0) { omFree(_rowKey); }
68  if (_numberOfColumnBlocks > 0) { omFree(_columnKey); }
69
70  _numberOfRowBlocks = lengthOfRowArray;
71  _numberOfColumnBlocks = lengthOfColumnArray;
72
73  /* allocate memory for new entries in _rowKey and _columnKey; */
74  _rowKey = (unsigned*)omAlloc(_numberOfRowBlocks*sizeof(unsigned));
75  _columnKey = (unsigned*)omAlloc(_numberOfColumnBlocks*sizeof(unsigned));
76
77  /* copying values from parameter arrays to private arrays */
78  for (int r = 0; r < _numberOfRowBlocks; r++)
79    _rowKey[r] = rowKey[r];
80  for (int c = 0; c < _numberOfColumnBlocks; c++)
81    _columnKey[c] = columnKey[c];
82}
83
84MinorKey::MinorKey(const int lengthOfRowArray,
85                   const unsigned int* const rowKey,
86                   const int lengthOfColumnArray,
87                   const unsigned int* const columnKey)
88{
89  _numberOfRowBlocks = lengthOfRowArray;
90  _numberOfColumnBlocks = lengthOfColumnArray;
91
92  /* allocate memory for new entries in _rowKey and _columnKey */
93  _rowKey = (unsigned*)omalloc(_numberOfRowBlocks*sizeof(unsigned));
94  _columnKey = (unsigned*)omalloc(_numberOfColumnBlocks*sizeof(unsigned));
95
96  /* copying values from parameter arrays to private arrays */
97  for (int r = 0; r < _numberOfRowBlocks; r++)
98    _rowKey[r] = rowKey[r];
99
100  for (int c = 0; c < _numberOfColumnBlocks; c++)
101    _columnKey[c] = columnKey[c];
102}
103
104MinorKey::~MinorKey()
105{
106  _numberOfRowBlocks = 0;
107  _numberOfColumnBlocks = 0;
108  omfree(_rowKey); _rowKey = NULL;
109  omfree(_columnKey); _columnKey = NULL;
110}
111
112//void MinorKey::print() const
113//{
114//  PrintS(this->toString().c_str());
115//}
116
117int MinorKey::getAbsoluteRowIndex(const int i) const
118{
119  /* This method is to return the absolute (0-based) index of the i-th
120     row encoded in \a this.
121     Example: bit-pattern of rows: "10010001101", i = 3:
122        This should yield the 0-based absolute index of the 3-rd bit
123        (counted from the right), i.e. 7. */
124
125  int matchedBits = -1; /* counter for matched bits;
126                           this needs to reach i, then we're done */
127  for (int block = 0; block < getNumberOfRowBlocks(); block ++)
128  {
129    /* start with lowest bits, i.e. in block No. 0 */
130    /* the bits in this block of 32 bits: */
131    unsigned int blockBits = getRowKey(block);
132    unsigned int shiftedBit = 1;
133    int exponent = 0;
134    /* The invariant "shiftedBit = 2^exponent" will hold throughout the
135       entire while loop. */
136    while (exponent < 32)
137    {
138      if (shiftedBit & blockBits) matchedBits++;
139      if (matchedBits == i) return exponent + (32 * block);
140      shiftedBit = shiftedBit << 1;
141      exponent++;
142    }
143  }
144  /* We should never reach this line of code. */
145  assume(false);
146  return -1;
147}
148
149int MinorKey::getAbsoluteColumnIndex(const int i) const
150{
151  /* This method is to return the absolute (0-based) index of the i-th
152     column encoded in \a this.
153     Example: bit-pattern of columns: "10010001101", i = 3:
154        This should yield the 0-based absolute index of the 3-rd bit
155        (counted from the right), i.e. 7. */
156
157  int matchedBits = -1; /* counter for matched bits; this needs to reach i,
158                           then we're done */
159  for (int block = 0; block < getNumberOfColumnBlocks(); block ++)
160  {
161    /* start with lowest bits, i.e. in block No. 0 */
162    /* the bits in this block of 32 bits: */
163    unsigned int blockBits = getColumnKey(block);
164    unsigned int shiftedBit = 1;
165    int exponent = 0;
166    /* The invariant "shiftedBit = 2^exponent" will hold throughout the
167       entire while loop. */
168    while (exponent < 32)
169    {
170      if (shiftedBit & blockBits) matchedBits++;
171      if (matchedBits == i) return exponent + (32 * block);
172      shiftedBit = shiftedBit << 1;
173      exponent++;
174    }
175  }
176  /* We should never reach this line of code. */
177  assume(false);
178  return -1;
179}
180
181void MinorKey::getAbsoluteRowIndices(int* const target) const
182{
183  int i = 0; /* index for filling the target array */
184  for (int block = 0; block < getNumberOfRowBlocks(); block ++)
185  {
186    /* start with lowest bits, i.e. in block No. 0 */
187    /* the bits in this block of 32 bits: */
188    unsigned int blockBits = getRowKey(block);
189    unsigned int shiftedBit = 1;
190    int exponent = 0;
191    /* The invariant "shiftedBit = 2^exponent" will hold throughout the
192       entire while loop. */
193    while (exponent < 32)
194    {
195      if (shiftedBit & blockBits) target[i++] = exponent + (32 * block);
196      shiftedBit = shiftedBit << 1;
197      exponent++;
198    }
199  }
200}
201
202void MinorKey::getAbsoluteColumnIndices(int* const target) const
203{
204  int i = 0; /* index for filling the target array */
205  for (int block = 0; block < getNumberOfColumnBlocks(); block ++)
206  {
207    /* start with lowest bits, i.e. in block No. 0 */
208    /* the bits in this block of 32 bits: */
209    unsigned int blockBits = getColumnKey(block);
210    unsigned int shiftedBit = 1;
211    int exponent = 0;
212    /* The invariant "shiftedBit = 2^exponent" will hold throughout the
213       entire while loop. */
214    while (exponent < 32)
215    {
216      if (shiftedBit & blockBits) target[i++] = exponent + (32 * block);
217      shiftedBit = shiftedBit << 1;
218      exponent++;
219    }
220  }
221}
222
223int MinorKey::getRelativeRowIndex(const int i) const
224{
225  /* This method is to return the relative (0-based) index of the row
226     with absolute index \c i.
227     Example: bit-pattern of rows: "10010001101", i = 7:
228        This should yield the 0-based relative index of the bit
229        corresponding to row no. 7, i.e. 3. */
230
231  int matchedBits = -1; /* counter for matched bits; this is going to
232                           contain our return value */
233  for (int block = 0; block < getNumberOfRowBlocks(); block ++)
234  {
235    /* start with lowest bits, i.e. in block No. 0 */
236    /* the bits in this block of 32 bits: */
237    unsigned int blockBits = getRowKey(block);
238    unsigned int shiftedBit = 1;
239    int exponent = 0;
240    /* The invariant "shiftedBit = 2^exponent" will hold throughout the
241       entire while loop. */
242    while (exponent < 32)
243    {
244      if (shiftedBit & blockBits) matchedBits++;
245      if (exponent + (32 * block) == i) return matchedBits;
246      shiftedBit = shiftedBit << 1;
247      exponent++;
248    }
249  }
250  /* We should never reach this line of code. */
251  assume(false);
252  return -1;
253}
254
255int MinorKey::getRelativeColumnIndex(const int i) const
256{
257  /* This method is to return the relative (0-based) index
258     of the column with absolute index \c i.
259     Example: bit-pattern of columns: "10010001101", i = 7:
260        This should yield the 0-based relative index of the bit
261        corresponding to column no. 7, i.e. 3. */
262
263  int matchedBits = -1; /* counter for matched bits; this is going
264                           to contain our return value */
265  for (int block = 0; block < getNumberOfColumnBlocks(); block ++)
266  {
267    /* start with lowest bits, i.e. in block No. 0 */
268    /* the bits in this block of 32 bits: */
269    unsigned int blockBits = getColumnKey(block);
270    unsigned int shiftedBit = 1;
271    int exponent = 0;
272    /* The invariant "shiftedBit = 2^exponent" will hold
273       throughout the entire while loop. */
274    while (exponent < 32)
275    {
276      if (shiftedBit & blockBits) matchedBits++;
277      if (exponent + (32 * block) == i) return matchedBits;
278      shiftedBit = shiftedBit << 1;
279      exponent++;
280    }
281  }
282  /* We should never reach this line of code. */
283  assume(false);
284  return -1;
285}
286
287unsigned int MinorKey::getRowKey(const int blockIndex) const
288{
289  return _rowKey[blockIndex];
290}
291
292unsigned int MinorKey::getColumnKey(const int blockIndex) const
293{
294  return _columnKey[blockIndex];
295}
296
297int MinorKey::getNumberOfRowBlocks() const
298{
299  return _numberOfRowBlocks;
300}
301
302int MinorKey::getNumberOfColumnBlocks() const
303{
304  return _numberOfColumnBlocks;
305}
306
307#ifndef SING_NDEBUG
308int MinorKey::getSetBits(const int a) const
309{
310  int b = 0;
311  if (a == 1)
312  { /* rows */
313    for (int i = 0; i < _numberOfRowBlocks; i++)
314    {
315      unsigned int m = _rowKey[i];
316      unsigned int k = 1;
317      for (int j = 0; j < 32; j++)
318      {
319        /* k = 2^j */
320        if (m & k) b++;
321        k = k << 1;
322      }
323    }
324  }
325  else
326  { /* columns */
327    for (int i = 0; i < _numberOfColumnBlocks; i++)
328    {
329      unsigned int m = _columnKey[i];
330      unsigned int k = 1;
331      for (int j = 0; j < 32; j++)
332      {
333        /* k = 2^j */
334        if (m & k) b++;
335        k = k << 1;
336      }
337    }
338  }
339  return b;
340}
341#endif
342
343MinorKey MinorKey::getSubMinorKey (const int absoluteEraseRowIndex,
344                                   const int absoluteEraseColumnIndex) const
345{
346  int rowBlock = absoluteEraseRowIndex / 32;
347  int exponent = absoluteEraseRowIndex % 32;
348  unsigned int newRowBits = getRowKey(rowBlock) - (1 << exponent);
349  int highestRowBlock = getNumberOfRowBlocks() - 1;
350  /* highestRowBlock will finally contain the highest block index with
351     non-zero bit pattern */
352  if ((newRowBits == 0) && (rowBlock == highestRowBlock))
353  {
354    /* we have thus nullified the highest block;
355       we can now forget about the highest block... */
356    highestRowBlock -= 1;
357    while (getRowKey(highestRowBlock) == 0) /* ...and maybe even some more
358                                               zero-blocks */
359      highestRowBlock -= 1;
360  }
361  /* highestRowBlock now contains the highest row block index with non-zero
362     bit pattern */
363
364  int columnBlock = absoluteEraseColumnIndex / 32;
365  exponent = absoluteEraseColumnIndex % 32;
366  unsigned int newColumnBits = getColumnKey(columnBlock) - (1 << exponent);
367  int highestColumnBlock = getNumberOfColumnBlocks() - 1;
368  /* highestColumnBlock will finally contain the highest block index with
369     non-zero bit pattern */
370  if ((newColumnBits == 0) && (columnBlock == highestColumnBlock))
371  {
372    /* we have thus nullified the highest block;
373       we can now forget about the highest block... */
374    highestColumnBlock -= 1;
375    while (getColumnKey(highestColumnBlock) == 0) /* ...and maybe even some
376                                                     more zero-blocks */
377      highestColumnBlock -= 1;
378  }
379  /* highestColumnBlock now contains the highest column block index with
380     non-zero bit pattern */
381
382  MinorKey result(highestRowBlock + 1, _rowKey, highestColumnBlock + 1,
383                  _columnKey);
384  /* This is just a copy with maybe some leading bit blocks omitted. We still
385     need to re-define the row block at index 'rowBlock' and the column block
386     at index 'columnBlock': */
387  if ((newRowBits != 0) || (rowBlock < getNumberOfRowBlocks() - 1))
388    result.setRowKey(rowBlock, newRowBits);
389  if ((newColumnBits != 0) || (columnBlock < getNumberOfColumnBlocks() - 1))
390    result.setColumnKey(columnBlock, newColumnBits);
391
392  #ifndef SING_NDEBUG
393  /* let's check that the number of selected rows and columns are equal;
394     (this check is only performed in the debug version) */
395  assume(result.getSetBits(1) == result.getSetBits(2));
396  #endif
397
398  return result;
399}
400
401void MinorKey::setRowKey (const int blockIndex, const unsigned int rowKey)
402{
403    _rowKey[blockIndex] = rowKey;
404}
405
406void MinorKey::setColumnKey (const int blockIndex,
407                             const unsigned int columnKey)
408{
409    _columnKey[blockIndex] = columnKey;
410}
411
412int MinorKey::compare (const MinorKey& that) const
413{
414  /* compare by rowKeys first; in case of equality, use columnKeys */
415  if (this->getNumberOfRowBlocks() < that.getNumberOfRowBlocks())
416    return -1;
417  if (this->getNumberOfRowBlocks() > that.getNumberOfRowBlocks())
418    return 1;
419  /* Here, numbers of rows are equal. */
420  for (int r = this->getNumberOfRowBlocks() - 1; r >= 0; r--)
421  {
422    if (this->getRowKey(r) < that.getRowKey(r)) return -1;
423    if (this->getRowKey(r) > that.getRowKey(r)) return 1;
424  }
425  /* Here, this and that encode ecaxtly the same sets of rows.
426     Now, we take a look at the columns. */
427  if (this->getNumberOfColumnBlocks() < that.getNumberOfColumnBlocks())
428    return -1;
429  if (this->getNumberOfColumnBlocks() > that.getNumberOfColumnBlocks())
430    return 1;
431  /* Here, numbers of columns are equal. */
432  for (int c = this->getNumberOfColumnBlocks() - 1; c >= 0; c--)
433  {
434    if (this->getColumnKey(c) < that.getColumnKey(c)) return -1;
435    if (this->getColumnKey(c) > that.getColumnKey(c)) return 1;
436  }
437  /* Here, this and that encode exactly the same sets of rows and columns. */
438  return 0;
439}
440
441/* just to make the compiler happy;
442   this method should never be called */
443bool MinorKey::operator==(const MinorKey& mk) const
444{
445  assume(false);
446  return this->compare(mk) == 0;
447}
448
449/* just to make the compiler happy;
450   this method should never be called */
451bool MinorKey::operator<(const MinorKey& mk) const
452{
453  assume(false);
454  return this->compare(mk) == -1;
455}
456
457void MinorKey::selectFirstRows (const int k, const MinorKey& mk)
458{
459  int hitBits = 0;      /* the number of bits we have hit; in the end, this
460                           has to be equal to k, the dimension of the minor */
461  int blockIndex = -1;  /* the index of the current int in mk */
462  unsigned int highestInt = 0;  /* the new highest block of this MinorKey */
463  /* We determine which ints of mk we can copy. Their indices will be
464     0, 1, ..., blockIndex - 1. And highestInt is going to capture the highest
465     int (which may be only a portion of the corresponding int in mk.
466     We loop until hitBits = k: */
467  while (hitBits < k)
468  {
469    blockIndex++;
470    highestInt = 0;
471    unsigned int currentInt = mk.getRowKey(blockIndex);
472    unsigned int shiftedBit = 1;
473    int exponent = 0;
474    /* invariant in the loop: shiftedBit = 2^exponent */
475    while (exponent < 32 && hitBits < k)
476    {
477      if (shiftedBit & currentInt)
478      {
479        highestInt += shiftedBit;
480        hitBits++;
481      }
482      shiftedBit = shiftedBit << 1;
483      exponent++;
484    }
485  }
486  /* free old memory */
487  omfree(_rowKey);
488  _rowKey = NULL;
489  _numberOfRowBlocks = blockIndex + 1;
490  /* allocate memory for new entries in _rowKey; */
491  _rowKey = (unsigned*)omAlloc(_numberOfRowBlocks*sizeof(unsigned));
492  /* copying values from mk to this MinorKey */
493  for (int r = 0; r < blockIndex; r++)
494    _rowKey[r] = mk.getRowKey(r);
495  _rowKey[blockIndex] = highestInt;
496}
497
498void MinorKey::selectFirstColumns (const int k, const MinorKey& mk)
499{
500  int hitBits = 0;      /* the number of bits we have hit; in the end, this
501                           has to be equal to k, the dimension of the minor */
502  int blockIndex = -1;  /* the index of the current int in mk */
503  unsigned int highestInt = 0;  /* the new highest block of this MinorKey */
504  /* We determine which ints of mk we can copy. Their indices will be
505     0, 1, ..., blockIndex - 1. And highestInt is going to capture the highest
506     int (which may be only a portion of the corresponding int in mk.
507     We loop until hitBits = k: */
508  while (hitBits < k)
509  {
510    blockIndex++;
511    highestInt = 0;
512    unsigned int currentInt = mk.getColumnKey(blockIndex);
513    unsigned int shiftedBit = 1;
514    int exponent = 0;
515    /* invariant in the loop: shiftedBit = 2^exponent */
516    while (exponent < 32 && hitBits < k)
517    {
518      if (shiftedBit & currentInt)
519      {
520        highestInt += shiftedBit;
521        hitBits++;
522      }
523      shiftedBit = shiftedBit << 1;
524      exponent++;
525    }
526  }
527  /* free old memory */
528  omfree(_columnKey); _columnKey = NULL;
529  _numberOfColumnBlocks = blockIndex + 1;
530  /* allocate memory for new entries in _columnKey; */
531  _columnKey = (unsigned*)omAlloc(_numberOfColumnBlocks*sizeof(unsigned));
532  /*  copying values from mk to this MinorKey */
533  for (int c = 0; c < blockIndex; c++)
534    _columnKey[c] = mk.getColumnKey(c);
535  _columnKey[blockIndex] = highestInt;
536}
537
538bool MinorKey::selectNextRows (const int k, const MinorKey& mk)
539{
540  /* We need to compute the set of k rows which must all be contained in mk.
541     AND: This set must be the least possible of this kind which is larger
542          than the currently encoded set of rows. (Here, '<' is w.r.t. to the
543          natural ordering on multi-indices.
544     Example: mk encodes the rows according to the bit pattern 11010111,
545              k = 3, this MinorKey encodes 10010100. Then, the method must
546              shift the set of rows in this MinorKey to 11000001 (, and
547              return true). */
548
549  /* The next two variables will finally name a row which is
550     (1) currently not yet among the rows in this MinorKey, but
551     (2) among the rows in mk, and
552     (3) which is "higher" than the lowest row in this MinorKey, and
553     (4) which is the lowest possible choice such that (1) - (3) hold.
554     If we should not be able to find such a row, then there is no next
555     subset of rows. In this case, the method will return false; otherwise
556     always true. */
557  int newBitBlockIndex = 0;        /* the block index of the bit */
558  unsigned int newBitToBeSet = 0;  /* the bit as 2^e, where 0 <= e <= 31 */
559
560  /* number of ints (representing rows) in this MinorKey: */
561  int blockCount = this->getNumberOfRowBlocks();
562  /* for iterating along the blocks of mk: */
563  int mkBlockIndex = mk.getNumberOfRowBlocks();
564
565  int hitBits = 0;    /* the number of bits we have hit */
566  int bitCounter = 0; /* for storing the number of bits hit before a
567                         specific moment; see below */
568  while (hitBits < k)
569  {
570    mkBlockIndex--;
571    unsigned int currentInt = mk.getRowKey(mkBlockIndex);
572    unsigned int shiftedBit = 1 << 31; /* initially, this equals 2^31, i.e.
573                                          the highest bit */
574    while (hitBits < k && shiftedBit > 0)
575    {
576      if ((blockCount - 1 >= mkBlockIndex) &&
577        (shiftedBit & this->getRowKey(mkBlockIndex))) hitBits++;
578      else if (shiftedBit & currentInt)
579      {
580        newBitToBeSet = shiftedBit;
581        newBitBlockIndex = mkBlockIndex;
582        bitCounter = hitBits; /* So, whenever we set newBitToBeSet, we want
583                                 to remember the momentary number of hit
584                                 bits. This will later be needed; see below. */
585      }
586      shiftedBit = shiftedBit >> 1;
587    }
588  }
589  if (newBitToBeSet == 0)
590  {
591    return false;
592  }
593  else
594  {
595    /* Note that the following must hold when reaching this line of code:
596       (1) The row with bit newBitToBeSet in this->getRowKey(newBitBlockIndex)
597           is currently not among the rows in this MinorKey, but
598       (2) it is among the rows in mk, and
599       (3) it is higher than the lowest row in this MinorKey, and
600       (4) it is the lowest possible choice such that (1) - (3) hold.
601       In the above example, we would reach this line with
602       newBitToBeSet == 2^6 and bitCounter == 1 (resulting from the bit 2^7).
603       */
604
605    if (blockCount - 1 < newBitBlockIndex)
606    { /* In this case, _rowKey is too small. */
607      /* free old memory */
608      omFree(_rowKey); _rowKey = NULL;
609      _numberOfRowBlocks = newBitBlockIndex + 1;
610      /* allocate memory for new entries in _rowKey; */
611      _rowKey = (unsigned*)omAlloc(_numberOfRowBlocks*sizeof(unsigned));
612      /* initializing entries to zero */
613        for (int r = 0; r < _numberOfRowBlocks; r++) _rowKey[r] = 0;
614    }
615    else
616    {
617      /* We need to delete all bits in _rowKey[newBitBlockIndex] that are
618         below newBitToBeSet: */
619      unsigned int anInt = this->getRowKey(newBitBlockIndex);
620      unsigned int deleteBit = newBitToBeSet >> 1; // in example: = 2^5
621      while (deleteBit > 0)
622      {
623        if (anInt & deleteBit) anInt -= deleteBit;
624        deleteBit = deleteBit >> 1;
625      };
626      _rowKey[newBitBlockIndex] = anInt;
627      /* ...and we delete all entries in _rowKey[i] for
628         0 <= i < newBitBlockIndex */
629      for (int i = 0; i < newBitBlockIndex; i++)
630        _rowKey[i] = 0;
631    }
632
633    /* We have now deleted all bits from _rowKey[...] below the bit
634       2^newBitToBeSet.
635       In the example we shall have at this point: _rowKey[...] = 10000000.
636       Now let's set the new bit: */
637    _rowKey[newBitBlockIndex] += newBitToBeSet;
638    /* in the example: _rowKey[newBitBlockIndex] = 11000000 */
639    bitCounter++; /* This is now the number of correct bits in _rowKey[...];
640                     i.e. in the example this will be equal to 2. */
641
642    /* Now we only need to fill _rowKey[...] with the lowest possible bits
643       until it consists of exactly k bits. (We know that we need to set
644       exactly (k - bitCounter) additional bits.) */
645    mkBlockIndex = -1;
646    while (bitCounter < k)
647    {
648      mkBlockIndex++;
649      unsigned int currentInt = mk.getRowKey(mkBlockIndex);
650      unsigned int shiftedBit = 1;
651      int exponent = 0;
652      /* invariant: shiftedBit = 2^exponent */
653      while (bitCounter < k && exponent < 32)
654      {
655        if (shiftedBit & currentInt)
656        {
657          _rowKey[mkBlockIndex] += shiftedBit;
658          bitCounter++;
659        };
660        shiftedBit = shiftedBit << 1;
661        exponent++;
662      }
663    };
664    /* in the example, we shall obtain _rowKey[...] = 11000001 */
665    return true;
666  }
667}
668
669bool MinorKey::selectNextColumns (const int k, const MinorKey& mk)
670{
671  /* We need to compute the set of k columns which must all be contained in mk.
672     AND: This set must be the least possible of this kind which is larger
673          than the currently encoded set of columns. (Here, '<' is w.r.t. to
674          the natural ordering on multi-indices.
675     Example: mk encodes the columns according to the bit pattern 11010111,
676              k = 3, this MinorKey encodes 10010100. Then, the method must
677              shift the set of columns in this MinorKey to 11000001 (, and
678              return true). */
679
680  /* The next two variables will finally name a column which is
681     (1) currently not yet among the columns in this MinorKey, but
682     (2) among the columns in mk, and
683     (3) which is "higher" than the lowest column in this MinorKey, and
684     (4) which is the lowest possible choice such that (1) - (3) hold.
685     If we should not be able to find such a column, then there is no next
686     subset of columns. In this case, the method will return false; otherwise
687     always true. */
688  int newBitBlockIndex = 0;        /* the block index of the bit */
689  unsigned int newBitToBeSet = 0;  /* the bit as 2^e, where 0 <= e <= 31 */
690
691   /* number of ints (representing columns) in this MinorKey: */
692  int blockCount = this->getNumberOfColumnBlocks();
693  /* for iterating along the blocks of mk: */
694  int mkBlockIndex = mk.getNumberOfColumnBlocks();
695
696  int hitBits = 0;    /* the number of bits we have hit */
697  int bitCounter = 0; /* for storing the number of bits hit before a specific
698                         moment; see below */
699  while (hitBits < k)
700  {
701    mkBlockIndex--;
702    unsigned int currentInt = mk.getColumnKey(mkBlockIndex);
703    unsigned int shiftedBit = 1 << 31; /* initially, this equals 2^31, i.e.
704                                          the highest bit */
705    while (hitBits < k && shiftedBit > 0)
706    {
707      if ((blockCount - 1 >= mkBlockIndex) &&
708        (shiftedBit & this->getColumnKey(mkBlockIndex))) hitBits++;
709      else if (shiftedBit & currentInt)
710      {
711        newBitToBeSet = shiftedBit;
712        newBitBlockIndex = mkBlockIndex;
713        bitCounter = hitBits; /* So, whenever we set newBitToBeSet, we want to
714                                 remember the momentary number of hit bits.
715                                 This will later be needed; see below. */
716      }
717      shiftedBit = shiftedBit >> 1;
718    }
719  }
720  if (newBitToBeSet == 0)
721  {
722    return false;
723  }
724  else
725  {
726    /* Note that the following must hold when reaching this line of code:
727       (1) The column with bit newBitToBeSet in
728           this->getColumnKey(newBitBlockIndex) is currently not among the
729           columns in this MinorKey, but
730       (2) it is among the columns in mk, and
731       (3) it is higher than the lowest columns in this MinorKey, and
732       (4) it is the lowest possible choice such that (1) - (3) hold.
733       In the above example, we would reach this line with
734       newBitToBeSet == 2^6 and bitCounter == 1 (resulting from the bit 2^7).
735       */
736
737    if (blockCount - 1 < newBitBlockIndex)
738    { /* In this case, _columnKey is too small. */
739        /* free old memory */
740        omFree( _columnKey); _columnKey = NULL;
741        _numberOfColumnBlocks = newBitBlockIndex + 1;
742        /* allocate memory for new entries in _columnKey; */
743        _columnKey = (unsigned*)omAlloc(_numberOfColumnBlocks*sizeof(unsigned));
744        /* initializing entries to zero */
745        for (int c = 0; c < _numberOfColumnBlocks; c++) _columnKey[c] = 0;
746    }
747    else
748    {
749      /* We need to delete all bits in _columnKey[newBitBlockIndex] that are
750         below newBitToBeSet: */
751      unsigned int anInt = this->getColumnKey(newBitBlockIndex);
752      unsigned int deleteBit = newBitToBeSet >> 1; /* in example: = 2^5 */
753      while (deleteBit > 0)
754      {
755        if (anInt & deleteBit) anInt -= deleteBit;
756        deleteBit = deleteBit >> 1;
757      };
758      _columnKey[newBitBlockIndex] = anInt;
759      /* ...and we delete all entries in _columnKey[i] fo
760         0 <= i < newBitBlockIndex */
761      for (int i = 0; i < newBitBlockIndex; i++)
762        _columnKey[i] = 0;
763    }
764    /* We have now deleted all bits from _columnKey[...] below the bit
765       2^newBitToBeSet. In the example we shall have at this point:
766       _columnKey[...] = 10000000. Now let's set the new bit: */
767    _columnKey[newBitBlockIndex] += newBitToBeSet;
768    /* in the example: _columnKey[newBitBlockIndex] = 11000000 */
769    bitCounter++; /* This is now the number of correct bits in
770                     _columnKey[...]; i.e. in the example this will be equal
771                     to 2. */
772
773    /* Now we only need to fill _columnKey[...] with the lowest possible bits
774       until it consists of exactly k bits. (We know that we need to set
775       exactly (k - bitCounter) additional bits.) */
776    mkBlockIndex = -1;
777    while (bitCounter < k)
778    {
779      mkBlockIndex++;
780      unsigned int currentInt = mk.getColumnKey(mkBlockIndex);
781      unsigned int shiftedBit = 1;
782      int exponent = 0;
783      /* invariant: shiftedBit = 2^exponent */
784      while (bitCounter < k && exponent < 32)
785      {
786        if (shiftedBit & currentInt)
787        {
788          _columnKey[mkBlockIndex] += shiftedBit;
789          bitCounter++;
790        };
791        shiftedBit = shiftedBit << 1;
792        exponent++;
793      }
794    };
795    /* in the example, we shall obtain _columnKey[...] = 11000001 */
796    return true;
797  }
798}
799
800string MinorKey::toString() const
801{ return ""; }
802/*
803  string t;
804  string s = "(";
805  unsigned int z = 0;
806  for (int r = this->getNumberOfRowBlocks() - 1; r >= 0; r--)
807  {
808    t = "";
809    z = this->getRowKey(r);
810    while (z != 0)
811    {
812      if ((z % 2) != 0) t = "1" + t; else t = "0" + t;
813      z = z / 2;
814    }
815    if (r < this->getNumberOfRowBlocks() - 1)
816      t = string(32 - t.length(), '0') + t;
817    s += t;
818  }
819  s += ", ";
820  for (int c = this->getNumberOfColumnBlocks() - 1; c >= 0; c--)
821  {
822    t = "";
823    z = this->getColumnKey(c);
824    while (z != 0)
825    {
826      if ((z % 2) != 0) t = "1" + t; else t = "0" + t;
827      z = z / 2;
828    }
829    if (c < this->getNumberOfColumnBlocks() - 1)
830      t = string(32 - t.length(), '0') + t;
831    s += t;
832  }
833  s += ")";
834  return s;
835}
836*/
837
838THREAD_VAR int MinorValue::g_rankingStrategy = -1;
839
840int MinorValue::getWeight () const
841{
842  assume(false);  /* must be overridden in derived classes */
843  return 0;
844}
845
846/* just to make the compiler happy;
847   this method should never be called */
848bool MinorValue::operator==(const MinorValue& mv) const
849{
850  assume(false);
851  return (this == &mv);  /* compare addresses of both objects */
852}
853
854string MinorValue::toString () const
855{
856  assume(false);  /* must be overridden in derived classes */
857  return "";
858}
859
860/* just to make the compiler happy;
861   this method should never be called */
862bool MinorValue::operator<(const MinorValue& mv) const
863{
864  assume(false);
865  return (this < &mv);  /* compare addresses of both objects */
866}
867
868int MinorValue::getRetrievals() const
869{
870  return _retrievals;
871}
872
873void MinorValue::incrementRetrievals()
874{
875  _retrievals++;
876}
877
878int MinorValue::getPotentialRetrievals() const
879{
880  return _potentialRetrievals;
881}
882
883int MinorValue::getMultiplications() const
884{
885  return _multiplications;
886}
887
888int MinorValue::getAdditions() const
889{
890  return _additions;
891}
892
893int MinorValue::getAccumulatedMultiplications() const
894{
895  return _accumulatedMult;
896}
897
898int MinorValue::getAccumulatedAdditions() const
899{
900  return _accumulatedSum;
901}
902
903void MinorValue::print() const
904{
905  PrintS(this->toString().c_str());
906}
907
908
909void MinorValue::SetRankingStrategy (const int rankingStrategy)
910{
911  g_rankingStrategy = rankingStrategy;
912  //if (g_rankingStrategy == 6) : rand() is never used
913  //{
914  //  /* initialize the random generator with system time */
915  //  srand ( time(NULL) );
916  //}
917}
918
919int MinorValue::GetRankingStrategy()
920{
921  return g_rankingStrategy;
922}
923
924/* this is for generically accessing the rank measure regardless of
925   which strategy has been set */
926int MinorValue::getUtility () const
927{
928  switch (this->GetRankingStrategy())
929  {
930    case 1: return this->rankMeasure1();
931    case 2: return this->rankMeasure2();
932    case 3: return this->rankMeasure3();
933    case 4: return this->rankMeasure4();
934    case 5: return this->rankMeasure5();
935    default: return this->rankMeasure1();
936  }
937}
938
939/* here are some sensible caching strategies: */
940int MinorValue::rankMeasure1 () const
941{
942  /* number of actually performed multiplications */
943  return this->getMultiplications();
944}
945
946int MinorValue::rankMeasure2 () const
947{
948  /* accumulated number of performed multiplications, i.e. all including
949     nested multiplications */
950  return this->getAccumulatedMultiplications();
951}
952
953int MinorValue::rankMeasure3 () const
954{
955  /* number of performed multiplications, weighted with the ratio of
956     not yet performed retrievals over the maximal number of retrievals */
957  return this->getMultiplications()
958         * (this->getPotentialRetrievals()
959         - this->getRetrievals())
960         / this->getPotentialRetrievals();
961}
962
963int MinorValue::rankMeasure4 () const
964{
965  /* number of performed multiplications,
966     multiplied with the number of not yet performed retrievals */
967  return this->getMultiplications()
968         * (this->getPotentialRetrievals()
969         - this->getRetrievals());
970}
971
972int MinorValue::rankMeasure5 () const
973{
974  /* number of not yet performed retrievals;
975     tends to cache entries longer when they are going to be retrieved more
976     often in the future */
977  return this->getPotentialRetrievals() - this->getRetrievals();
978}
979
980int IntMinorValue::getWeight () const
981{
982  /* put measure for size of MinorValue here, i.e. number of monomials in
983     polynomial; so far, we use the accumulated number of multiplications
984     (i.e., including all nested ones) to simmulate the size of a polynomial */
985  return _accumulatedMult;
986}
987
988IntMinorValue::IntMinorValue (const int result, const int multiplications,
989                              const int additions,
990                              const int accumulatedMultiplications,
991                              const int accumulatedAdditions,
992                              const int retrievals,
993                              const int potentialRetrievals)
994{
995  _result = result;
996  _multiplications = multiplications;
997  _additions = additions;
998  _accumulatedMult = accumulatedMultiplications;
999  _accumulatedSum = accumulatedAdditions;
1000  _potentialRetrievals = potentialRetrievals;
1001  _retrievals = retrievals;
1002}
1003
1004IntMinorValue::IntMinorValue ()
1005{
1006  _result = -1;
1007  _multiplications = -1;
1008  _additions = -1;
1009  _accumulatedMult = -1;
1010  _accumulatedSum = -1;
1011  _potentialRetrievals = -1;
1012  _retrievals = -1;
1013}
1014
1015IntMinorValue::~IntMinorValue()
1016{
1017}
1018
1019int IntMinorValue::getResult() const
1020{
1021  return _result;
1022}
1023
1024string IntMinorValue::toString () const
1025{
1026  char h[10];
1027
1028  /* Let's see whether a cache has been used to compute this MinorValue: */
1029  bool cacheHasBeenUsed = true;
1030  if (this->getRetrievals() == -1) cacheHasBeenUsed = false;
1031
1032  sprintf(h, "%d", this->getResult());
1033  string s = h;
1034  s += " [retrievals: ";
1035  if (cacheHasBeenUsed) { sprintf(h, "%d", this->getRetrievals()); s += h; }
1036  else s += "/";
1037  s += " (of ";
1038  if (cacheHasBeenUsed)
1039  {
1040    sprintf(h, "%d", this->getPotentialRetrievals());
1041    s += h;
1042  }
1043  else s += "/";
1044  s += "), *: ";
1045  sprintf(h, "%d", this->getMultiplications()); s += h;
1046  s += " (accumulated: ";
1047  sprintf(h, "%d", this->getAccumulatedMultiplications()); s += h;
1048  s += "), +: ";
1049  sprintf(h, "%d", this->getAdditions()); s += h;
1050  s += " (accumulated: ";
1051  sprintf(h, "%d", this->getAccumulatedAdditions()); s += h;
1052  s += "), rank: ";
1053  if (cacheHasBeenUsed) { sprintf(h, "%d", this->getUtility()); s += h; }
1054  else s += "/";
1055  s += "]";
1056  return s;
1057}
1058
1059IntMinorValue::IntMinorValue (const IntMinorValue& mv)
1060{
1061  _result = mv.getResult();
1062  _retrievals = mv.getRetrievals();
1063  _potentialRetrievals = mv.getPotentialRetrievals();
1064  _multiplications = mv.getMultiplications();
1065  _additions = mv.getAdditions();
1066  _accumulatedMult = mv.getAccumulatedMultiplications();
1067  _accumulatedSum = mv.getAccumulatedAdditions();
1068}
1069
1070PolyMinorValue::PolyMinorValue (const poly result, const int multiplications,
1071                                const int additions,
1072                                const int accumulatedMultiplications,
1073                                const int accumulatedAdditions,
1074                                const int retrievals,
1075                                const int potentialRetrievals)
1076{
1077  _result = pCopy(result);
1078  _multiplications = multiplications;
1079  _additions = additions;
1080  _accumulatedMult = accumulatedMultiplications;
1081  _accumulatedSum = accumulatedAdditions;
1082  _potentialRetrievals = potentialRetrievals;
1083  _retrievals = retrievals;
1084}
1085
1086PolyMinorValue::PolyMinorValue ()
1087{
1088  _result = NULL;
1089  _multiplications = -1;
1090  _additions = -1;
1091  _accumulatedMult = -1;
1092  _accumulatedSum = -1;
1093  _potentialRetrievals = -1;
1094  _retrievals = -1;
1095}
1096
1097PolyMinorValue::~PolyMinorValue()
1098{
1099  p_Delete(&_result, currRing);
1100}
1101
1102poly PolyMinorValue::getResult() const
1103{
1104  return _result;
1105}
1106
1107int PolyMinorValue::getWeight () const
1108{
1109  /* put measure for size of PolyMinorValue here, e.g. the number of monomials
1110     in the cached polynomial */
1111  return pLength(_result); // the number of monomials in the polynomial
1112}
1113
1114string PolyMinorValue::toString () const
1115{
1116  char h[20];
1117
1118  /* Let's see whether a cache has been used to compute this MinorValue: */
1119  bool cacheHasBeenUsed = true;
1120  if (this->getRetrievals() == -1) cacheHasBeenUsed = false;
1121
1122  string s = pString(_result);
1123  s += " [retrievals: ";
1124  if (cacheHasBeenUsed) { sprintf(h, "%d", this->getRetrievals()); s += h; }
1125  else s += "/";
1126  s += " (of ";
1127  if (cacheHasBeenUsed)
1128  {
1129    sprintf(h, "%d", this->getPotentialRetrievals());
1130    s += h;
1131  }
1132  else s += "/";
1133  s += "), *: ";
1134  sprintf(h, "%d", this->getMultiplications()); s += h;
1135  s += " (accumulated: ";
1136  sprintf(h, "%d", this->getAccumulatedMultiplications()); s += h;
1137  s += "), +: ";
1138  sprintf(h, "%d", this->getAdditions()); s += h;
1139  s += " (accumulated: ";
1140  sprintf(h, "%d", this->getAccumulatedAdditions()); s += h;
1141  s += "), rank: ";
1142  if (cacheHasBeenUsed) { sprintf(h, "%d", this->getUtility()); s += h; }
1143  else s += "/";
1144  s += "]";
1145  return s;
1146}
1147
1148PolyMinorValue::PolyMinorValue (const PolyMinorValue& mv)
1149{
1150  _result = pCopy(mv.getResult());
1151  _retrievals = mv.getRetrievals();
1152  _potentialRetrievals = mv.getPotentialRetrievals();
1153  _multiplications = mv.getMultiplications();
1154  _additions = mv.getAdditions();
1155  _accumulatedMult = mv.getAccumulatedMultiplications();
1156  _accumulatedSum = mv.getAccumulatedAdditions();
1157}
1158
1159void PolyMinorValue::operator= (const PolyMinorValue& mv)
1160{
1161  if (_result != mv.getResult()) pDelete(&_result);
1162  _result = pCopy(mv.getResult());
1163  _retrievals = mv.getRetrievals();
1164  _potentialRetrievals = mv.getPotentialRetrievals();
1165  _multiplications = mv.getMultiplications();
1166  _additions = mv.getAdditions();
1167  _accumulatedMult = mv.getAccumulatedMultiplications();
1168  _accumulatedSum = mv.getAccumulatedAdditions();
1169}
Note: See TracBrowser for help on using the repository browser.