1 /* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2007 Bjarke Hammersholt Roune (www.broune.com)
3 Copyright (C) 2010 University of Aarhus
4 Contact Bjarke Hammersholt Roune for license information (www.broune.com)
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program. If not, see http://www.gnu.org/licenses/.
18 */
19 #include "stdinc.h"
20 #include "IdealTree.h"
21
22 #include "Ideal.h"
23 #include "Term.h"
24
25 namespace {
26 const size_t MaxLeafSize = 60;
27
rawStrictlyDivides(Ideal::const_iterator begin,Ideal::const_iterator end,const Exponent * term,size_t varCount)28 bool rawStrictlyDivides(Ideal::const_iterator begin,
29 Ideal::const_iterator end,
30 const Exponent* term,
31 size_t varCount) {
32 for (; begin != end; ++begin)
33 if (Term::strictlyDivides(*begin, term, varCount))
34 return true;
35 return false;
36 }
37 }
38
39 class IdealTree::Node {
40 public:
Node(Ideal::iterator begin,Ideal::iterator end,size_t varCount)41 Node(Ideal::iterator begin,
42 Ideal::iterator end,
43 size_t varCount):
44 _begin(begin), _end(end), _varCount(varCount) {
45 makeTree();
46 }
47
48 void makeTree();
49 bool strictlyContains(const Exponent* term) const;
getVarCount() const50 size_t getVarCount() const {return _varCount;}
51
52 private:
53 auto_ptr<Node> _lessOrEqual;
54 auto_ptr<Node> _greater;
55 Ideal::iterator _begin;
56 Ideal::iterator _end;
57 size_t _varCount;
58 size_t _var;
59 size_t _pivot;
60 };
61
makeTree()62 void IdealTree::Node::makeTree() {
63 if ((size_t)distance(_begin, _end) <= MaxLeafSize)
64 return;
65
66 Term lcm(_varCount);
67 Term gcd(*_begin, _varCount);
68 for (Ideal::const_iterator it = _begin; it != _end; ++it) {
69 lcm.lcm(lcm, *it);
70 gcd.gcd(gcd, *it);
71 }
72
73 while (true) {
74 size_t maxVar = 0;
75 for (size_t var = 1; var < _varCount; ++var)
76 if (lcm[var] - gcd[var] > lcm[maxVar] - gcd[maxVar])
77 maxVar = var;
78
79 // TODO: could this ever happen?
80 if (lcm[maxVar] == 0)
81 break; // we are not making any progress anyway.
82
83 ASSERT(lcm[maxVar] >= 1);
84 _var = maxVar;
85
86 // It's significant that we are rounding down to ensure that
87 // neither of _lessOrEqual and _greater becomes empty.
88 _pivot = (lcm[maxVar] + gcd[maxVar]) >> 1; // Note: x >> 1 == x / 2
89
90 Ideal::iterator left = _begin;
91 Ideal::iterator right = _end - 1;
92
93 // put those strictly divisible by _var^_pivot to right and the
94 // rest to the left.
95 while (left != right) {
96 // Find left-most element that should go right.
97 while ((*left)[_var] <= _pivot && left != right)
98 ++left;
99
100 // Find right-most element that should go left.
101 while ((*right)[_var] > _pivot && left != right)
102 --right;
103
104 // Swap the two found elements so that they both go into the
105 // right place.
106 using std::swap;
107 swap(*left, *right);
108 }
109 ASSERT((*(_end - 1))[_var] > _pivot);
110 ASSERT((*_begin)[_var] <= _pivot);
111
112 // Make middle the first element that went right, so that the two
113 // ranges become [_begin, middle) and [middle, _end).
114 ASSERT((*right)[_var] > _pivot)
115 Ideal::iterator middle = right;
116 while ((*middle)[_var] > _pivot) {
117 ASSERT(middle != _begin);
118 --middle;
119 }
120 ++middle;
121 ASSERT(_begin < middle && middle <= _end);
122 ASSERT((*middle)[_var] > _pivot);
123 ASSERT((*(middle - 1))[_var] <= _pivot);
124
125 _lessOrEqual.reset(new Node(_begin, middle, _varCount));
126 _greater.reset(new Node(middle, _end, _varCount));
127
128 _lessOrEqual->makeTree();
129 _greater->makeTree();
130 return;
131 }
132 }
133
strictlyContains(const Exponent * term) const134 bool IdealTree::Node::strictlyContains(const Exponent* term) const {
135 if (_lessOrEqual.get() != 0) {
136 ASSERT(_greater.get() != 0);
137 bool returnValue =
138 _lessOrEqual->strictlyContains(term) ||
139 _greater->strictlyContains(term);
140 ASSERT(returnValue == rawStrictlyDivides(_begin, _end, term, _varCount));
141 return returnValue;
142 } else {
143 ASSERT(_greater.get() == 0);
144 return rawStrictlyDivides(_begin, _end, term, _varCount);
145 }
146 }
147
IdealTree(const Ideal & ideal)148 IdealTree::IdealTree(const Ideal& ideal) {
149 // not using initialization for this to avoid depending on order of
150 // initialization of members.
151 _storage.reset(new Ideal(ideal));
152 _root.reset
153 (new Node(_storage->begin(), _storage->end(), ideal.getVarCount()));
154 }
155
~IdealTree()156 IdealTree::~IdealTree() {
157 // Destructor defined so auto_ptr<T> in the header does not need
158 // definition of T.
159 }
160
strictlyContains(const Exponent * term) const161 bool IdealTree::strictlyContains(const Exponent* term) const {
162 ASSERT(_root.get() != 0);
163 return _root->strictlyContains(term);
164 }
165
getVarCount() const166 size_t IdealTree::getVarCount() const {
167 ASSERT(_root.get() != 0);
168 return _root->getVarCount();
169 }
170