Frobby  0.9.0
RawSquareFreeIdeal.cpp
Go to the documentation of this file.
1 /* Frobby: Software for monomial ideal computations.
2  Copyright (C) 2010 University of Aarhus
3  Contact Bjarke Hammersholt Roune for license information (www.broune.com)
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see http://www.gnu.org/licenses/.
17 */
18 #include "stdinc.h"
19 #include "RawSquareFreeIdeal.h"
20 
21 #include "Arena.h"
22 #include "Ideal.h"
23 #include "RawSquareFreeTerm.h"
24 #include "BigIdeal.h"
25 #include <limits>
26 #include <algorithm>
27 #include <sstream>
28 #include <cstring>
29 
31 namespace Ops = SquareFreeTermOps;
32 
33 namespace {
41  template<class Iter, class Pred>
42  inline Iter RsfPartition(Iter begin, Iter end, Pred pred, size_t varCount) {
43  // The invariant of the loop is that pred(x) is true if x precedes
44  // begin and pred(x) is false if x is at or after end.
45  while (true) {
46  if (begin == end)
47  return begin;
48  while (pred(*begin))
49  if (++begin == end)
50  return begin;
51  // now pred(*begin) is true and begin < end.
52  if (begin == --end)
53  return begin;
54  while (!pred(*end))
55  if (begin == --end)
56  return begin;
57  // now pred(*end) is false and begin < end.
58 
59  // This swap is the reason that std::partition doesn't work
60  Ops::swap(*begin, *end, varCount);
61  ++begin;
62  }
63  }
64 
70  const size_t wordCount) {
71  for (RSFIdeal::iterator it = begin; it != end;) {
72  for (RSFIdeal::const_iterator div = begin; div != end; ++div) {
73  if (Ops::divides(*div, *div + wordCount, *it) && div != it) {
74  --end;
75  Ops::assign(*it, *it + wordCount, *end);
76  goto next;
77  }
78  }
79  ++it;
80  next:;
81  }
82  return end;
83  }
84 }
85 
86 RSFIdeal* RSFIdeal::construct(void* buffer, size_t varCount) {
87  RSFIdeal* p = static_cast<RSFIdeal*>(buffer);
88  p->_varCount = varCount;
89  p->_wordsPerTerm = Ops::getWordCount(varCount);
90  p->_genCount = 0;
91  p->_memoryEnd = p->_memory;
92  ASSERT(p->isValid());
93  return p;
94 }
95 
96 RSFIdeal* RSFIdeal::construct(void* buffer, const Ideal& ideal) {
97  RSFIdeal* p = construct(buffer, ideal.getVarCount());
98  p->insert(ideal);
99  ASSERT(p->isValid());
100  return p;
101 }
102 
103 RSFIdeal* RSFIdeal::construct(void* buffer, const RawSquareFreeIdeal& ideal) {
104  RSFIdeal* p = construct(buffer, ideal.getVarCount());
105  p->insert(ideal);
106  ASSERT(p->isValid());
107  return p;
108 }
109 
110 size_t RSFIdeal::getBytesOfMemoryFor(size_t varCount, size_t generatorCount) {
111  // This calculation is tricky because there are many overflows that
112  // can occur. If most cases if an overflow occurs or nearly occurs
113  // then the amount of memory needed could not be allocated on the
114  // system. In this case 0 is returned to signal the error. Note that
115  // x / y rounded up is (x - 1) / y + 1 for x, y > 0. The
116  // multiplication a * b does not overflow if a <= MAX /
117  // b. Otherwise, there may not be an overflow, but a * b will still
118  // be too much to reasonably allocate.
119 
120  size_t bytesForStruct = sizeof(RSFIdeal) - sizeof(Word);
121  if (generatorCount == 0)
122  return bytesForStruct;
123 
124  // Compute bytes per generator taking into account memory alignment.
125  size_t bytesPerGenUnaligned = varCount == 0 ? 1 : (varCount - 1) / 8 + 1;
126  size_t wordsPerGen = (bytesPerGenUnaligned - 1) / sizeof(Word) + 1;
127  if (wordsPerGen > numeric_limits<size_t>::max() / sizeof(Word))
128  return 0;
129  size_t bytesPerGen = wordsPerGen * sizeof(Word);
130 
131  // Compute bytes in all.
132  if (bytesPerGen > numeric_limits<size_t>::max() / generatorCount)
133  return 0;
134  size_t bytesForGens = bytesPerGen * generatorCount;
135  if (bytesForGens > numeric_limits<size_t>::max() - bytesForStruct)
136  return 0;
137  return bytesForStruct + bytesForGens;
138 }
139 
140 void RSFIdeal::setToTransposeOf(const RawSquareFreeIdeal& ideal, Word* eraseVars) {
141  if (this == &ideal) {
142  transpose(eraseVars);
143  return;
144  }
145 
146  const size_t idealVarCount = ideal.getVarCount();
147  const size_t idealGenCount = ideal.getGeneratorCount();
148 
149  _varCount = idealGenCount;
151  _genCount = 0;
153 
154  for (size_t var = 0; var < idealVarCount; ++var) {
155  if (eraseVars != 0 && Ops::getExponent(eraseVars, var))
156  continue;
157  insertIdentity();
158  Word* newTransposedGen = back();
159  for (size_t gen = 0; gen < idealGenCount; ++gen) {
160  const bool value = Ops::getExponent(ideal.getGenerator(gen), var);
161  Ops::setExponent(newTransposedGen, gen, value);
162  }
163  }
164 
165  ASSERT(isValid());
166 }
167 
168 void RSFIdeal::transpose(Word* eraseVars) {
169  const size_t varCount = getVarCount();
170  const size_t genCount = getGeneratorCount();
171  const size_t bytes = RSFIdeal::getBytesOfMemoryFor(varCount, genCount);
172  Arena& arena = Arena::getArena();
173  void* buffer = arena.alloc(bytes);
174  RSFIdeal* copy = RSFIdeal::construct(buffer, *this);
175  setToTransposeOf(*copy, eraseVars);
176  arena.freeTop(buffer);
177 }
178 
179 void RSFIdeal::print(FILE* file) const {
180  ostringstream out;
181  print(out);
182  fputs(out.str().c_str(), file);
183 }
184 
185 void RSFIdeal::print(ostream& out) const {
186  const size_t varCount = getVarCount();
187  out << "//------------ Ideal (Square Free):\n";
188  for (size_t gen = 0; gen < getGeneratorCount(); ++gen) {
189  for (size_t var = 0; var < varCount; ++var)
190  out << Ops::getExponent(getGenerator(gen), var);
191  out << '\n';
192  }
193  out << "------------\\\\\n";
194 }
195 
196 size_t RSFIdeal::insert(const Ideal& ideal) {
197  ASSERT(getVarCount() == ideal.getVarCount());
198 
199  size_t gen = 0;
200  for (; gen < ideal.getGeneratorCount(); ++gen) {
201  if (!Ops::encodeTerm(_memoryEnd, ideal[gen], getVarCount()))
202  break;
203  ++_genCount;
205  }
206  ASSERT(isValid());
207  return gen;
208 }
209 
210 size_t RSFIdeal::insert(const BigIdeal& ideal) {
211  ASSERT(getVarCount() == ideal.getVarCount());
212 
213  size_t gen = 0;
214  for (; gen < ideal.getGeneratorCount(); ++gen) {
215  if (!Ops::encodeTerm(_memoryEnd, ideal[gen], getVarCount()))
216  break;
217  ++_genCount;
219  }
220  ASSERT(isValid());
221  return gen;
222 }
223 
225  const_iterator stop = ideal.end();
226  for (const_iterator it = ideal.begin(); it != stop; ++it)
227  insert(*it);
228  ASSERT(isValid());
229 }
230 
231 bool RSFIdeal::insert(const std::vector<std::string>& term) {
232  ASSERT(term.size() == getVarCount());
233 
234  if (!Ops::encodeTerm(_memoryEnd, term, getVarCount()))
235  return false;
236  ++_genCount;
238  ASSERT(isValid());
239  return true;
240 }
241 
243  iterator newEnd = ::minimize(begin(), end(), getWordsPerTerm());
244  _genCount = newEnd - begin();
245  _memoryEnd = *newEnd;
246  ASSERT(isValid());
247 }
248 
249 void RSFIdeal::colon(const Word* by) {
250  const size_t wordCount = getWordsPerTerm();
251  const iterator stop = end();
252  for (iterator it = begin(); it != stop; ++it)
253  Ops::colonInPlace(*it, *it + wordCount, by);
254  ASSERT(isValid());
255 }
256 
257 void RSFIdeal::colon(size_t var) {
258  const iterator stop = end();
259  for (iterator it = begin(); it != stop; ++it)
260  Ops::setExponent(*it, var, false);
261  ASSERT(isValid());
262 }
263 
264 void RSFIdeal::compact(const Word* remove) {
265  const size_t oldVarCount = getVarCount();
266  const iterator oldBegin = begin();
267  const iterator oldStop = end();
268 
269  // Compact each term without moving that term.
270  size_t varCompact = 0;
271  for (size_t var = 0; var < oldVarCount; ++var) {
272  if (Ops::getExponent(remove, var) != 0)
273  continue;
274  for (iterator oldIt = oldBegin; oldIt != oldStop; ++oldIt)
275  Ops::setExponent(*oldIt, varCompact, Ops::getExponent(*oldIt, var));
276  ++varCompact;
277  }
278  // varCompact is now the number of variables in the compacted ideal.
279 
280  // The last word in each term must have zeroes at those positions
281  // that are past the number of actual variables. So we need to go through and
282  const size_t bitOffset = Ops::getBitOffset(varCompact);
283  const size_t wordOffset = Ops::getWordOffset(varCompact);
284  if (bitOffset != 0) {
285  const Word mask = (((Word)1) << bitOffset) - 1;
286  for (iterator oldIt = oldBegin; oldIt != oldStop; ++oldIt)
287  *(*oldIt + wordOffset) &= mask;
288  }
289 
290  // Copy the new compacted terms to remove the space between them. We
291  // couldn't do that before because those spaces contained exponents
292  // that we had not extracted yet.
293  const size_t newWordCount = Ops::getWordCount(varCompact);
294  iterator newIt(_memory, newWordCount);
295  for (iterator oldIt = oldBegin; oldIt != oldStop; ++oldIt, ++newIt)
296  Ops::assign(*newIt, (*newIt) + newWordCount, *oldIt);
297 
298  _varCount = varCompact;
299  _wordsPerTerm = newWordCount;
300  _memoryEnd = *newIt;
301  ASSERT(isValid());
302 }
303 
304 void RSFIdeal::getLcmOfNonMultiples(Word* lcm, size_t var) const {
305  ASSERT(var < getVarCount());
306 
307  const Word* const lcmEnd = lcm + getWordsPerTerm();
308  Ops::setToIdentity(lcm, lcmEnd);
309 
310  const const_iterator stop = end();
311  for (const_iterator it = begin(); it != stop; ++it)
312  if (Ops::getExponent(*it, var) == 0)
313  Ops::lcmInPlace(lcm, lcmEnd, *it);
314  ASSERT(isValid());
315 }
316 
317 static inline void countVarDividesBlockUpTo15(const Word* it,
318  size_t genCount,
319  const size_t wordsPerTerm,
320  size_t* counts) {
321  // mask has the bit pattern 0001 0001 ... 0001
322  const Word mask = ~((Word)0u) / 15u;
323 
324  Word a0, a1, a2, a3;
325  if ((genCount & 1) == 1) {
326  const Word a = *it;
327  a0 = a & mask;
328  a1 = (a >> 1) & mask;
329  a2 = (a >> 2) & mask;
330  a3 = (a >> 3) & mask;
331  it += wordsPerTerm;
332  } else
333  a0 = a1 = a2 = a3 = 0;
334 
335  genCount >>= 1;
336  for (size_t i = 0; i < genCount; ++i) {
337  const Word a = *it;
338  it += wordsPerTerm;
339  const Word aa = *it;
340  it += wordsPerTerm;
341 
342  a0 += a & mask;
343  a1 += (a >> 1) & mask;
344  a2 += (a >> 2) & mask;
345  a3 += (a >> 3) & mask;
346 
347  a0 += aa & mask;
348  a1 += (aa >> 1) & mask;
349  a2 += (aa >> 2) & mask;
350  a3 += (aa >> 3) & mask;
351  }
352 
353  for (size_t i = 0; i < BitsPerWord / 4; ++i) {
354  *(counts + 0) += a0 & 0xF;
355  *(counts + 1) += a1 & 0xF;
356  *(counts + 2) += a2 & 0xF;
357  *(counts + 3) += a3 & 0xF;
358  a0 >>= 4;
359  a1 >>= 4;
360  a2 >>= 4;
361  a3 >>= 4;
362  counts += 4;
363  }
364 }
365 
366 void RSFIdeal::getVarDividesCounts(vector<size_t>& divCounts) const {
367  const size_t varCount = getVarCount();
368  const size_t wordCount = getWordsPerTerm();
369  // We reserve BitsPerWord extra space. Otherwise we would have to
370  // make sure not to index past the end of the vector of counts when
371  // dealing with variables in the unused part of the last word.
372  divCounts.reserve(getVarCount() + BitsPerWord);
373  divCounts.resize(getVarCount());
374  size_t* divCountsBasePtr = &(divCounts.front());
375  size_t* divCountsEnd = divCountsBasePtr + BitsPerWord * wordCount;
376  memset(divCountsBasePtr, 0, sizeof(size_t) * varCount);
377 
378  size_t generatorsToGo = getGeneratorCount();
379  const_iterator blockBegin = begin();
380  while (generatorsToGo > 0) {
381  const size_t blockSize = generatorsToGo >= 15 ? 15 : generatorsToGo;
382 
383  size_t* counts = divCountsBasePtr;
384  const Word* genOffset = *blockBegin;
385  for (; counts != divCountsEnd; counts += BitsPerWord, ++genOffset)
386  countVarDividesBlockUpTo15(genOffset, blockSize, wordCount, counts);
387 
388  generatorsToGo -= blockSize;
389  blockBegin = blockBegin + blockSize;
390  }
391 }
392 
393 size_t RSFIdeal::getMultiple(size_t var) const {
394  ASSERT(var < getVarCount());
395 
396  const const_iterator stop = end();
397  const const_iterator start = begin();
398  for (const_iterator it = start; it != stop; ++it)
399  if (Ops::getExponent(*it, var) == 1)
400  return it - start;
401  return getGeneratorCount();
402 }
403 
404 size_t RSFIdeal::getNonMultiple(size_t var) const {
405  ASSERT(var < getVarCount());
406 
407  const const_iterator stop = end();
408  const const_iterator start = begin();
409  for (const_iterator it = start; it != stop; ++it)
410  if (Ops::getExponent(*it, var) == 0)
411  return it - start;
412  return getGeneratorCount();
413 }
414 
416  const_iterator it = begin();
417  const const_iterator stop = end();
418  if (it == stop)
419  return 0;
420 
421  const size_t varCount = getVarCount();
422  size_t maxSupp = Ops::getSizeOfSupport(*it, varCount);
423  const_iterator maxSuppIt = it;
424 
425  for (++it; it != stop; ++it) {
426  const size_t supp = Ops::getSizeOfSupport(*it, varCount);
427  if (maxSupp < supp) {
428  maxSupp = supp;
429  maxSuppIt = it;
430  }
431  }
432  return maxSuppIt - begin();
433 }
434 
436  const_iterator it = begin();
437  const const_iterator stop = end();
438  if (it == stop)
439  return 0;
440 
441  const size_t varCount = getVarCount();
442  size_t minSupp = Ops::getSizeOfSupport(*it, varCount);
443  const_iterator minSuppIt = it;
444 
445  for (++it; it != stop; ++it) {
446  const size_t supp = Ops::getSizeOfSupport(*it, varCount);
447  if (minSupp > supp) {
448  minSupp = supp;
449  minSuppIt = it;
450  }
451  }
452  return minSuppIt - begin();
453 }
454 
455 void RSFIdeal::getGcdOfMultiples(Word* gcd, size_t var) const {
456  ASSERT(var < getVarCount());
457 
458  const Word* const gcdEnd = gcd + getWordsPerTerm();
460 
461  const const_iterator stop = end();
462  for (const_iterator it = begin(); it != stop; ++it)
463  if (Ops::getExponent(*it, var) == 1)
464  Ops::gcdInPlace(gcd, gcdEnd, *it);
465 }
466 
467 void RSFIdeal::getGcdOfMultiples(Word* gcd, const Word* div) const {
468  const size_t varCount = getVarCount();
469  const size_t wordCount = getWordsPerTerm();
470  const Word* const gcdEnd = gcd + wordCount;
471  const Word* divEnd = div + wordCount;
472  Ops::setToAllVarProd(gcd, varCount);
473 
474  const const_iterator stop = end();
475  for (const_iterator it = begin(); it != stop; ++it)
476  if (Ops::divides(div, divEnd, *it))
477  Ops::gcdInPlace(gcd, gcdEnd, *it);
478 }
479 
480 void RSFIdeal::removeGenerator(size_t gen) {
481  Word* term = getGenerator(gen);
482  Word* last = _memoryEnd - getWordsPerTerm();
483  if (term != last)
484  Ops::assign(term, term + getWordsPerTerm(), last);
485  --_genCount;
487  ASSERT(isValid());
488 }
489 
491  const RawSquareFreeIdeal& ideal) {
492  ASSERT(getVarCount() == ideal.getVarCount());
493 
494  const Word* termEnd = term + getWordsPerTerm();
495  const_iterator stop = ideal.end();
496  for (const_iterator it = ideal.begin(); it != stop; ++it)
497  if (!Ops::divides(term, termEnd, *it))
498  insert(*it);
499  ASSERT(isValid());
500 }
501 
503  const RawSquareFreeIdeal& ideal) {
504  ASSERT(getVarCount() == ideal.getVarCount());
505 
506  const_iterator stop = ideal.end();
507  for (const_iterator it = ideal.begin(); it != stop; ++it)
508  if (Ops::getExponent(*it, var) == 0)
509  insert(*it);
510  ASSERT(isValid());
511 }
512 
514  const size_t wordCount = getWordsPerTerm();
515  for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
516  if (!Ops::isRelativelyPrime(term, term + wordCount, getGenerator(gen)))
517  return gen;
518  return getGeneratorCount();
519 }
520 
522  const size_t wordCount = getWordsPerTerm();
523  for (size_t offset = 0; offset < wordCount; ++offset) {
524  Word once = 0;
525  Word twice = 0;
526  for (size_t gen = 0; gen < _genCount; ++gen) {
527  Word word = getGenerator(gen)[offset];
528  twice |= once & word;
529  once |= word;
530  }
531  const Word onceOnly = once & ~twice;
532  if (onceOnly != 0) {
533  for (size_t gen = 0; ; ++gen) {
534  ASSERT(gen < _genCount);
535  Word word = getGenerator(gen)[offset];
536  if (word & onceOnly)
537  return gen;
538  }
539  ASSERT(false);
540  }
541  }
542  return getGeneratorCount();
543 }
544 
545 bool RSFIdeal::hasFullSupport(const Word* ignore) const {
546  ASSERT(ignore != 0);
547  const size_t wordCount = getWordsPerTerm();
548  const Word allOnes = ~((Word)0);
549 
550  const Word* firstGenOffset = _memory;
551  const Word* endGenOffset = _memoryEnd;
552  size_t varsLeft = getVarCount();
553  while (true) {
554  const Word* gen = firstGenOffset;
555  Word support = *ignore;
556  while (gen != endGenOffset) {
557  support |= *gen;
558  gen += wordCount;
559  }
560 
561  if (varsLeft > BitsPerWord) {
562  if (support != allOnes)
563  return false;
564  } else {
565  if (varsLeft == BitsPerWord)
566  return support == allOnes;
567  const Word fullSupportWord = (((Word)1) << varsLeft) - 1;
568  return support == fullSupportWord;
569  }
570 
571  varsLeft -= BitsPerWord;
572  ++ignore;
573  ++firstGenOffset;
574  ++endGenOffset;
575  }
576  return true;
577 }
578 
580  size_t wordCount = getWordsPerTerm();
581  for (size_t i = 0; i < _genCount; ++i)
582  for (size_t div = 0; div < _genCount; ++div)
583  if (div != i &&
584  Ops::divides(getGenerator(div), getGenerator(div) + wordCount,
585  getGenerator(i)))
586  return false;
587  return true;
588 }
589 
590 void RSFIdeal::swap(size_t a, size_t b) {
591  ASSERT(a < getGeneratorCount());
592  ASSERT(b < getGeneratorCount());
594  ASSERT(isValid());
595 }
596 
597 bool RSFIdeal::operator==(const RawSquareFreeIdeal& ideal) const {
598  if (getVarCount() != ideal.getVarCount())
599  return false;
600  if (getGeneratorCount() != ideal.getGeneratorCount())
601  return false;
602 
603  const size_t varCount = getVarCount();
604  const_iterator stop = end();
605  const_iterator it = begin();
606  const_iterator it2 = ideal.begin();
607  for (; it != stop; ++it, ++it2)
608  if (!Ops::equals(*it, *it2, varCount))
609  return false;
610  return true;
611 }
612 
613 namespace {
614  struct CmpForSortLexAscending : std::binary_function<size_t, size_t, bool> {
615  bool operator()(size_t a, size_t b) const {
616  return Ops::lexLess(ideal->getGenerator(a), ideal->getGenerator(b),
617  ideal->getVarCount());
618  }
619  RawSquareFreeIdeal* ideal;
620  };
621 }
622 
624  vector<size_t> sorted(getGeneratorCount());
625  for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
626  sorted[gen] = gen;
627  {
628  CmpForSortLexAscending cmp;
629  cmp.ideal = this;
630  std::sort(sorted.begin(), sorted.end(), cmp);
631  }
632 
633  RawSquareFreeIdeal* clone =
635  for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
636  clone->insert(getGenerator(gen));
637  for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
638  Ops::assign(getGenerator(gen), clone->getGenerator(sorted[gen]),
639  getVarCount());
641  ASSERT(isValid());
642 }
643 
644 void RSFIdeal::insert(const Word* term) {
646  ++_genCount;
648  ASSERT(isValid());
649 }
650 
652  const iterator stop = end();
653  const size_t varCount = getVarCount();
654  for (iterator it = begin(); it != stop; ++it)
655  Ops::invert(*it, varCount);
656 }
657 
660  ++_genCount;
662  ASSERT(isValid());
663 }
664 
665 namespace {
666  struct ColonReminimizeTermHelper {
667  bool operator()(const Word* term) {
668  return !Ops::isRelativelyPrime(colon, colonEnd, term);
669  }
670  const Word* colon;
671  const Word* colonEnd;
672  };
673 }
674 
676  ASSERT(by != 0);
677  const size_t varCount = getVarCount();
678  const size_t wordCount = getWordsPerTerm();
679  const iterator start = begin();
680  iterator stop = end();
681 
682  ColonReminimizeTermHelper colonIsRelativelyPrime;
683  colonIsRelativelyPrime.colon = by;
684  colonIsRelativelyPrime.colonEnd = by + wordCount;
685  iterator middle = RsfPartition(start, stop, colonIsRelativelyPrime, varCount);
686 
687  if (middle == start) {
688  ASSERT(isValid());
689  return; // colon is relatively prime to all generators
690  }
691  for (iterator it = start; it != middle; ++it)
692  Ops::colonInPlace(*it, *it + wordCount, by);
693 
694  // var is not relatively prime to [start, middle) and is relatively
695  // prime to [middle, end).
696 
697  iterator newMiddle = ::minimize(start, middle, getWordsPerTerm());
698 
699  iterator newEnd = newMiddle;
700  for (iterator it = middle; it != stop; ++it) {
701  for (const_iterator div = start; div != newMiddle; ++div)
702  if (Ops::divides(*div, *div + wordCount, *it))
703  goto next;
704  Ops::assign(*newEnd, *newEnd + wordCount, *it);
705  ++newEnd;
706  next:;
707  }
708 
709  _memoryEnd = *newEnd;
710  _genCount = newEnd - start;
711  ASSERT(isValid());
712 }
713 
714 namespace {
715  struct ColonReminimizeVarHelper {
716  bool operator()(const Word* term) {
717  return Ops::getExponent(term, var) != 0;
718  }
719  size_t var;
720  };
721 }
722 
723 void RSFIdeal::colonReminimize(size_t var) {
724  ASSERT(var < getVarCount());
725  const size_t varCount = getVarCount();
726  const size_t wordCount = getWordsPerTerm();
727  const iterator start = begin();
728  iterator stop = end();
729 
730  ColonReminimizeVarHelper varDivides;
731  varDivides.var = var;
732  iterator middle = RsfPartition(start, stop, varDivides, varCount);
733 
734  if (middle == start) {
735  ASSERT(isValid());
736  return; // var divides no generators
737  }
738  for (iterator it = start; it != middle; ++it)
739  Ops::setExponent(*it, var, 0);
740  if (middle == stop) {
741  ASSERT(isValid());
742  return; // var divides all
743  }
744 
745  // var divided [start, middle) and did (does) not divide [middle,
746  // end). Both of these ranges are minimized on their own, and no
747  // element of [middle, end) divides an element of [start, middle).
748  for (iterator it = middle; it != stop;) {
749  for (const_iterator div = start; div != middle; ++div) {
750  if (Ops::divides(*div, *div + wordCount, *it)) {
751  --stop;
752  Ops::assign(*it, *it + wordCount, *stop);
753  --_genCount;
754  goto next;
755  }
756  }
757  ++it;
758  next:;
759  }
760  _memoryEnd = *stop;
761 
762  ASSERT(isValid());
763 }
764 
765 void RSFIdeal::getLcm(Word* lcm) const {
766  Word* lcmEnd = lcm + getWordsPerTerm();
767  Ops::setToIdentity(lcm, lcmEnd);
768  const const_iterator stop = end();
769  for (const_iterator it = begin(); it != stop; ++it)
770  Ops::lcmInPlace(lcm, lcmEnd, *it);
771  ASSERT(isValid());
772 }
773 
774 bool RSFIdeal::isValid() const {
775  const size_t varCount = getVarCount();
776  const size_t wordCount = getWordsPerTerm();
777  const size_t genCount = getGeneratorCount();
778  if (varCount != _varCount)
779  return false;
780  if (wordCount != _wordsPerTerm)
781  return false;
782  if (_genCount != genCount)
783  return false;
784 
785  if (wordCount != Ops::getWordCount(varCount))
786  return false;
787  if (_memoryEnd != _memory + wordCount * genCount)
788  return false;
789  if (_memoryEnd < _memory)
790  return false; // happens on overflow
791 
792  for (const Word* p = _memory; p != _memoryEnd; p += wordCount)
793  Ops::isValid(p, varCount);
794  return true;
795 }
796 
797 RSFIdeal* newRawSquareFreeIdeal(size_t varCount, size_t capacity) {
798  size_t byteCount = RSFIdeal::getBytesOfMemoryFor(varCount, capacity);
799  if (byteCount == 0)
800  throw bad_alloc();
801  void* buffer = new char[byteCount];
802  RSFIdeal* ideal = RSFIdeal::construct(buffer, varCount);
803  ASSERT(ideal->isValid());
804  return ideal;
805 }
806 
808  size_t byteCount = RSFIdeal::getBytesOfMemoryFor(ideal.getVarCount(),
809  ideal.getGeneratorCount());
810  if (byteCount == 0)
811  throw bad_alloc();
812  void* buffer = new char[byteCount];
813  RawSquareFreeIdeal* p = RSFIdeal::construct(buffer, ideal.getVarCount());
814  p->insert(ideal);
815  ASSERT(p->isValid());
816  return p;
817 }
818 
820  istringstream in(str);
821  vector<string> lines;
822  string line;
823 
824  while (getline(in, line))
825  if (line != "")
826  lines.push_back(line);
827 
828  const size_t varCount = lines.empty() ? 0 : lines.front().size();
829  RawSquareFreeIdeal* ideal = newRawSquareFreeIdeal(varCount, lines.size());
830  for (size_t gen = 0; gen < lines.size(); ++gen) {
831  ASSERT(lines[gen].size() == varCount);
832  Word* term = Ops::newTermParse(lines[gen].c_str());
833  ideal->insert(term);
834  Ops::deleteTerm(term);
835  }
836  ASSERT(ideal->isValid());
837  return ideal;
838 }
839 
841  delete[] reinterpret_cast<char*>(ideal);
842 }
void sortLexAscending()
Sorts the generators in ascending lex order.
RawSquareFreeIdeal RSFIdeal
size_t getMaxSupportGen() const
Returns the index of a generator with maximum support.
void lcmInPlace(Word *res, const Word *resEnd, const Word *a)
void getVarDividesCounts(vector< size_t > &counts) const
Sets counts[var] to the number of generators that var divides.
size_t getExclusiveVarGenerator()
Returns the index of a generator that is the only one to be divisible by some variable.
bool lexLess(const Word *a, const Word *b, size_t varCount)
bool equals(const Word *a, const Word *b, size_t varCount)
Returns true if a equals b.
const_iterator doesn't have all it needs to be a proper STL iterator.
size_t getMinSupportGen() const
Returns the index of a generator with minimum support.
void colonReminimize(const Word *colon)
Performs a colon and minimize.
bool hasFullSupport(const Word *ignore) const
Returns true if for every variable it either divides ignore or it divides some (not necessarily minim...
void invert(Word *a, size_t varCount)
Make 0 exponents 1 and make 1 exponents 0.
size_t getWordOffset(size_t var)
Word * getGenerator(size_t index)
Returns the generator at index.
#define ASSERT(X)
Definition: stdinc.h:85
void removeGenerator(size_t index)
Removes the generator at index.
Represents a monomial ideal with int exponents.
Definition: Ideal.h:27
size_t getBitOffset(size_t var)
void lcm(Word *res, const Word *resEnd, const Word *a, const Word *b)
void transpose(Word *eraseVars=0)
Equivalent to setToTransposeOf(this, eraseVars).
RSFIdeal * newRawSquareFreeIdeal(size_t varCount, size_t capacity)
Allocates object with enough memory for capacity generators in varCount variables.
size_t getGeneratorCount() const
Definition: BigIdeal.h:144
void setToIdentity(Word *res, const Word *resEnd)
size_t getNotRelativelyPrime(const Word *term)
Returns the index of the first generator that is not relatively prime with term.
A bit packed square free ideal placed in a pre-allocated buffer.
void setToTransposeOf(const RawSquareFreeIdeal &ideal, Word *eraseVars=0)
Resets this object to the transpose of ideal.
void colon(Word *res, const Word *resEnd, const Word *a, const Word *b)
void gcd(Word *res, const Word *resEnd, const Word *a, const Word *b)
void setExponent(Word *a, size_t var, bool value)
size_t getGeneratorCount() const
size_t getVarCount() const
Definition: Ideal.h:56
void compact(const Word *remove)
Removes the variables that divide remove.
iterator doesn't have all it needs to be a proper STL iterator.
void getLcmOfNonMultiples(Word *lcm, size_t var) const
Sets lcm to be the least common multple of those generators that var does not divide.
static Arena & getArena()
Returns an arena object that can be used for non-thread safe scratch memory after static objects have...
Definition: Arena.h:126
unsigned long Word
The native unsigned type for the CPU.
Definition: stdinc.h:92
static RawSquareFreeIdeal * construct(void *buffer, size_t varCount=0)
This is an arena allocator.
Definition: Arena.h:53
static const size_t BitsPerWord
Definition: stdinc.h:93
void * alloc(size_t size)
Returns a pointer to a buffer of size bytes.
Definition: Arena.h:180
static void countVarDividesBlockUpTo15(const Word *it, size_t genCount, const size_t wordsPerTerm, size_t *counts)
void swap(size_t a, size_t b)
bool encodeTerm(Word *encoded, const Exponent *term, const size_t varCount)
Assigns the RawSquareFreeTerm-encoded form of term to encoded and returns true if term is square free...
bool getExponent(const Word *a, size_t var)
returns true if var divides a and false otherwise.
size_t getWordsPerTerm() const
size_t getWordCount(size_t varCount)
void print(FILE *file) const
Print a debug-suitable representation of this object to file.
bool isValid(const Word *a, size_t varCount)
The unused bits at the end of the last word must be zero for the functions here to work correctly...
void swap01Exponents()
Change 0 exponents into 1 and vice versa.
void freeTop(void *ptr)
Frees the buffer pointed to by ptr.
Definition: Arena.h:209
bool operator==(const RawSquareFreeIdeal &ideal) const
Returns true if *this equals ideal.
RawSquareFreeIdeal * newRawSquareFreeIdealParse(const char *str)
Allocates and returns an ideal based on str.
void colonInPlace(Word *res, const Word *resEnd, const Word *b)
bool isRelativelyPrime(const Word *a, const Word *b, size_t varCount)
void setToAllVarProd(Word *res, size_t varCount)
Sets all exponents of res to 1.
size_t getSizeOfSupport(const Word *a, size_t varCount)
Word * newTermParse(const char *strParam)
Allocates and returns a term based on str.
size_t getNonMultiple(size_t var) const
Returns the index of the first generator that var does not divide or getGeneratorCount() if no such g...
bool isMinimallyGenerated() const
Returns true if no generator divides another.
void assign(Word *a, const Word *b, size_t varCount)
void getLcm(Word *lcm) const
Puts the least common multiple of the generators of the ideal into lcm.
size_t getVarCount() const
Definition: BigIdeal.h:148
bool isValid() const
Returns true if the internal invariants of ideal are satisfied.
size_t insert(const Ideal &ideal)
Inserts the generators of ideal from index 0 onward until reaching a non-squarefree generator or all ...
size_t getVarCount() const
size_t getMultiple(size_t var) const
Returns the index of the first generator that var divides or getGeneratorCount() if no such generator...
void insertNonMultiples(const Word *term, const RawSquareFreeIdeal &ideal)
Insert those generators of ideal that are not multiples of term.
void deleteRawSquareFreeIdeal(RSFIdeal *ideal)
static size_t getBytesOfMemoryFor(size_t varCount, size_t generatorCount)
Returns the number of bytes of memory necessary to contain an ideal with the given parameters...
void colon(const Word *by)
bool divides(const Word *a, const Word *aEnd, const Word *b)
Returns true if a divides b.
void getGcdOfMultiples(Word *gcd, size_t var) const
Sets gcd to be the greatest common denominator of those generators that are divisible by var...
void deleteTerm(Word *term)
Deletes term previously returned by newTerm().
void swap(hashtable< _Val, _Key, _HF, _Extract, _EqKey, _All > &__ht1, hashtable< _Val, _Key, _HF, _Extract, _EqKey, _All > &__ht2)
Definition: hashtable.h:740
void gcdInPlace(Word *res, const Word *resEnd, const Word *a)
size_t getGeneratorCount() const
Definition: Ideal.h:57