Frobby 0.9.5
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
31namespace Ops = SquareFreeTermOps;
32
33namespace {
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
86RSFIdeal* 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
96RSFIdeal* 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
103RSFIdeal* 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
110size_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
140void 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;
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
168void 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
179void RSFIdeal::print(FILE* file) const {
180 ostringstream out;
181 print(out);
182 fputs(out.str().c_str(), file);
183}
184
185void 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
196size_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
210size_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
231bool RSFIdeal::insert(const std::vector<std::string>& term) {
232 ASSERT(term.size() == getVarCount());
233
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
249void 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
257void 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
264void 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
304void 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
317static 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
366void 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
393size_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
404size_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
455void 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
467void 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
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
545bool 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
590void RSFIdeal::swap(size_t a, size_t b) {
594 ASSERT(isValid());
595}
596
597bool 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
613namespace {
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
644void 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
664
665namespace {
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
714namespace {
715 struct ColonReminimizeVarHelper {
716 bool operator()(const Word* term) {
717 return Ops::getExponent(term, var) != 0;
718 }
719 size_t var;
720 };
721}
722
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
765void 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
774bool 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
797RSFIdeal* 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];
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}
RSFIdeal * newRawSquareFreeIdeal(size_t varCount, size_t capacity)
Allocates object with enough memory for capacity generators in varCount variables.
RawSquareFreeIdeal RSFIdeal
void deleteRawSquareFreeIdeal(RSFIdeal *ideal)
RawSquareFreeIdeal * newRawSquareFreeIdealParse(const char *str)
Allocates and returns an ideal based on str.
static void countVarDividesBlockUpTo15(const Word *it, size_t genCount, const size_t wordsPerTerm, size_t *counts)
RawSquareFreeIdeal * newRawSquareFreeIdeal(size_t varCount, size_t capacity)
Allocates object with enough memory for capacity generators in varCount variables.
void deleteRawSquareFreeIdeal(RawSquareFreeIdeal *ideal)
Deallocates memory returned by newRawSquareFreeIdeal().
This is an arena allocator.
Definition Arena.h:53
void * alloc(size_t size)
Returns a pointer to a buffer of size bytes.
Definition Arena.h:180
void freeTop(void *ptr)
Frees the buffer pointed to by ptr.
Definition Arena.h:209
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
size_t getVarCount() const
Definition BigIdeal.h:148
size_t getGeneratorCount() const
Definition BigIdeal.h:144
Represents a monomial ideal with int exponents.
Definition Ideal.h:27
size_t getGeneratorCount() const
Definition Ideal.h:57
size_t getVarCount() const
Definition Ideal.h:56
const_iterator doesn't have all it needs to be a proper STL iterator.
iterator doesn't have all it needs to be a proper STL iterator.
A bit packed square free ideal placed in a pre-allocated buffer.
void sortLexAscending()
Sorts the generators in ascending lex order.
static RawSquareFreeIdeal * construct(void *buffer, size_t varCount=0)
size_t getNotRelativelyPrime(const Word *term)
Returns the index of the first generator that is not relatively prime with term.
bool isValid() const
Returns true if the internal invariants of ideal are satisfied.
bool hasFullSupport(const Word *ignore) const
Returns true if for every variable it either divides ignore or it divides some (not necessarily minim...
void getGcdOfMultiples(Word *gcd, size_t var) const
Sets gcd to be the greatest common denominator of those generators that are divisible by var.
bool operator==(const RawSquareFreeIdeal &ideal) const
Returns true if *this equals ideal.
void colon(const Word *by)
size_t getExclusiveVarGenerator()
Returns the index of a generator that is the only one to be divisible by some variable.
void getLcm(Word *lcm) const
Puts the least common multiple of the generators of the ideal into lcm.
size_t getMinSupportGen() const
Returns the index of a generator with minimum support.
void getVarDividesCounts(vector< size_t > &counts) const
Sets counts[var] to the number of generators that var divides.
void compact(const Word *remove)
Removes the variables that divide remove.
void removeGenerator(size_t index)
Removes the generator at index.
void print(FILE *file) const
Print a debug-suitable representation of this object to file.
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.
size_t getMaxSupportGen() const
Returns the index of a generator with maximum support.
size_t getVarCount() const
Word * getGenerator(size_t index)
Returns the generator at index.
void colonReminimize(const Word *colon)
Performs a colon and minimize.
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...
size_t getMultiple(size_t var) const
Returns the index of the first generator that var divides or getGeneratorCount() if no such generator...
size_t getGeneratorCount() const
size_t insert(const Ideal &ideal)
Inserts the generators of ideal from index 0 onward until reaching a non-squarefree generator or all ...
void transpose(Word *eraseVars=0)
Equivalent to setToTransposeOf(this, eraseVars).
void swap01Exponents()
Change 0 exponents into 1 and vice versa.
void swap(size_t a, size_t b)
void setToTransposeOf(const RawSquareFreeIdeal &ideal, Word *eraseVars=0)
Resets this object to the transpose of ideal.
bool isMinimallyGenerated() const
Returns true if no generator divides another.
void getLcmOfNonMultiples(Word *lcm, size_t var) const
Sets lcm to be the least common multple of those generators that var does not divide.
void insertNonMultiples(const Word *term, const RawSquareFreeIdeal &ideal)
Insert those generators of ideal that are not multiples of term.
size_t getWordsPerTerm() const
size_t getWordOffset(size_t var)
void setExponent(Word *a, size_t var, bool value)
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.
size_t getWordCount(size_t varCount)
void colon(Word *res, const Word *resEnd, const Word *a, const Word *b)
size_t getBitOffset(size_t var)
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...
Word * newTermParse(const char *strParam)
Allocates and returns a term based on str.
void invert(Word *a, size_t varCount)
Make 0 exponents 1 and make 1 exponents 0.
bool lexLess(const Word *a, const Word *b, size_t varCount)
size_t getSizeOfSupport(const Word *a, size_t varCount)
void colonInPlace(Word *res, const Word *resEnd, const Word *b)
void lcmInPlace(Word *res, const Word *resEnd, const Word *a)
bool getExponent(const Word *a, size_t var)
returns true if var divides a and false otherwise.
void setToAllVarProd(Word *res, size_t varCount)
Sets all exponents of res to 1.
bool equals(const Word *a, const Word *b, size_t varCount)
Returns true if a equals b.
bool divides(const Word *a, const Word *aEnd, const Word *b)
Returns true if a divides b.
void swap(Word *a, Word *b, size_t varCount)
void deleteTerm(Word *term)
Deletes term previously returned by newTerm().
void gcdInPlace(Word *res, const Word *resEnd, const Word *a)
void assign(Word *a, const Word *b, size_t varCount)
bool isRelativelyPrime(const Word *a, const Word *b, size_t varCount)
void setToIdentity(Word *res, const Word *resEnd)
This header file includes common definitions and is included as the first line of code in every imple...
static const size_t BitsPerWord
Definition stdinc.h:94
unsigned long Word
The native unsigned type for the CPU.
Definition stdinc.h:93
#define ASSERT(X)
Definition stdinc.h:86