source: git/Singular/Minor.cc @ d2ea299

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