1 | #include "mod2.h" |
---|
2 | #include <stdio.h> |
---|
3 | #include <string> |
---|
4 | #include <kernel/matpol.h> |
---|
5 | #include <kernel/intvec.h> |
---|
6 | #include <Singular/lists.h> |
---|
7 | #include <Cone.h> |
---|
8 | #include <Fan.h> |
---|
9 | |
---|
10 | #ifdef HAVE_FANS |
---|
11 | |
---|
12 | Fan::Fan(const intvec* maxRays, |
---|
13 | const intvec* facetNs, |
---|
14 | const intvec* linSpace) |
---|
15 | { |
---|
16 | m_maxRays = NULL; |
---|
17 | if (maxRays != NULL) m_maxRays = ivCopy(maxRays); |
---|
18 | m_facetNs = NULL; |
---|
19 | if (facetNs != NULL) m_facetNs = ivCopy(facetNs); |
---|
20 | m_linSpace = NULL; |
---|
21 | if (linSpace != NULL) m_linSpace = ivCopy(linSpace); |
---|
22 | else |
---|
23 | { |
---|
24 | /* trivial lineality space: we use a matrix with a single zero column */ |
---|
25 | if (m_maxRays != NULL) m_linSpace = new intvec(m_maxRays->rows(), 1, 0); |
---|
26 | else m_linSpace = new intvec(m_facetNs->rows(), 1, 0); |
---|
27 | } |
---|
28 | m_maxCones = (lists)omAllocBin(slists_bin); |
---|
29 | m_maxCones->Init(0); |
---|
30 | m_dim = -1; |
---|
31 | m_complete = -1; |
---|
32 | m_simplicial = -1; |
---|
33 | m_pure = -1; |
---|
34 | m_adjacency = (lists)omAllocBin(slists_bin); |
---|
35 | m_adjacency->Init(0); |
---|
36 | } |
---|
37 | |
---|
38 | intvec* Fan::getMaxRays() const { return m_maxRays; } |
---|
39 | intvec* Fan::getFacetNs() const { return m_facetNs; } |
---|
40 | intvec* Fan::getLinSpace() const { return m_linSpace; } |
---|
41 | lists Fan::getMaxCones() const { return m_maxCones; } |
---|
42 | int Fan::getDim() const { return m_dim; } |
---|
43 | int Fan::getComplete() const { return m_complete; } |
---|
44 | int Fan::getSimplicial() const { return m_simplicial; } |
---|
45 | int Fan::getPure() const { return m_pure; } |
---|
46 | lists Fan::getRawAdj() const { return m_adjacency; } |
---|
47 | |
---|
48 | void Fan::setDim (const int dim) |
---|
49 | { |
---|
50 | m_dim = dim; |
---|
51 | if (m_dim < 0) m_dim = -1; |
---|
52 | } |
---|
53 | |
---|
54 | void Fan::setComplete (const int complete) |
---|
55 | { |
---|
56 | m_complete = complete; |
---|
57 | if ((m_complete != 0) && (m_complete != 1)) m_complete = -1; |
---|
58 | } |
---|
59 | |
---|
60 | void Fan::setSimplicial (const int simplicial) |
---|
61 | { |
---|
62 | m_simplicial = simplicial; |
---|
63 | if ((m_simplicial != 0) && (m_simplicial != 1)) m_simplicial = -1; |
---|
64 | } |
---|
65 | |
---|
66 | void Fan::setPure (const int pure) |
---|
67 | { |
---|
68 | m_pure = pure; |
---|
69 | if ((m_pure != 0) && (m_pure != 1)) m_pure = -1; |
---|
70 | } |
---|
71 | |
---|
72 | int Fan::getAmbientDim() const |
---|
73 | { |
---|
74 | if (m_linSpace == NULL) return -1; |
---|
75 | else return m_linSpace->rows(); |
---|
76 | } |
---|
77 | |
---|
78 | Fan::Fan(const Fan& f) |
---|
79 | { /* this is deep copy */ |
---|
80 | m_maxRays = NULL; |
---|
81 | if (f.getMaxRays() != NULL) m_maxRays = ivCopy(f.getMaxRays()); |
---|
82 | m_facetNs = NULL; |
---|
83 | if (f.getFacetNs() != NULL) m_facetNs = ivCopy(f.getFacetNs()); |
---|
84 | m_linSpace = ivCopy(f.getLinSpace()); |
---|
85 | m_maxCones = lCopy(f.getMaxCones()); |
---|
86 | m_dim = f.getDim(); |
---|
87 | m_complete = f.getComplete(); |
---|
88 | m_simplicial = f.getSimplicial(); |
---|
89 | m_pure = f.getPure(); |
---|
90 | m_adjacency = lCopy(f.getRawAdj()); |
---|
91 | } |
---|
92 | |
---|
93 | Fan::Fan() |
---|
94 | { |
---|
95 | m_maxRays = NULL; |
---|
96 | m_facetNs = NULL; |
---|
97 | m_linSpace = NULL; |
---|
98 | m_maxCones = (lists)omAllocBin(slists_bin); |
---|
99 | m_maxCones->Init(0); |
---|
100 | m_dim = -1; |
---|
101 | m_complete = -1; |
---|
102 | m_simplicial = -1; |
---|
103 | m_pure = -1; |
---|
104 | m_adjacency = (lists)omAllocBin(slists_bin); |
---|
105 | m_adjacency->Init(0); |
---|
106 | } |
---|
107 | |
---|
108 | Fan::~Fan() |
---|
109 | { |
---|
110 | if (m_maxRays != NULL) delete m_maxRays; |
---|
111 | if (m_facetNs != NULL) delete m_facetNs; |
---|
112 | if (m_linSpace != NULL) delete m_linSpace; |
---|
113 | if (m_maxCones->nr + 1 > 0) |
---|
114 | omFreeSize((ADDRESS)m_maxCones->m, |
---|
115 | (m_maxCones->nr + 1)*sizeof(sleftv)); |
---|
116 | omFreeBin((ADDRESS)m_maxCones, slists_bin); |
---|
117 | if (m_adjacency->nr + 1 > 0) |
---|
118 | omFreeSize((ADDRESS)m_adjacency->m, |
---|
119 | (m_adjacency->nr + 1)*sizeof(sleftv)); |
---|
120 | omFreeBin((ADDRESS)m_adjacency, slists_bin); |
---|
121 | } |
---|
122 | |
---|
123 | char* Fan::toString() const |
---|
124 | { |
---|
125 | char h[50]; |
---|
126 | std::string s; |
---|
127 | if (m_linSpace == NULL) |
---|
128 | s = "// fan (no properties set)"; |
---|
129 | else |
---|
130 | { |
---|
131 | s = "// fan:\n"; |
---|
132 | int i = getAmbientDim(); |
---|
133 | sprintf(h, "// ambient dimension: %d\n", i); s += h; |
---|
134 | i = 0; if (m_maxRays != NULL) i = m_maxRays->cols(); |
---|
135 | sprintf(h, "// %d maximal rays\n", i); s += h; |
---|
136 | i = 0; if (m_facetNs != NULL) i = m_facetNs->cols(); |
---|
137 | sprintf(h, "// %d facet normals\n", i); s += h; |
---|
138 | i = m_maxCones->nr + 1; |
---|
139 | sprintf(h, "// %d maximal cones\n", i); s += h; |
---|
140 | s += "// dimension: " ; |
---|
141 | if (m_dim == -1) s += "?\n"; |
---|
142 | else { sprintf(h, "%d\n", m_dim); s += h; } |
---|
143 | s += "// complete: " ; |
---|
144 | if (m_complete == -1) s += "?\n"; |
---|
145 | else if (m_complete == 0) s += "no\n"; |
---|
146 | else s += "yes\n"; |
---|
147 | s += "// simplicial: " ; |
---|
148 | if (m_simplicial == -1) s += "?\n"; |
---|
149 | else if (m_simplicial == 0) s += "no\n"; |
---|
150 | else s += "yes\n"; |
---|
151 | s += "// pure: " ; |
---|
152 | if (m_pure == -1) s += "?"; |
---|
153 | else if (m_pure == 0) s += "no"; |
---|
154 | else s += "yes"; |
---|
155 | } |
---|
156 | return omStrDup(s.c_str()); /* any caller must clean up this string */ |
---|
157 | } |
---|
158 | |
---|
159 | void Fan::print() const |
---|
160 | { |
---|
161 | char* s = this->toString(); |
---|
162 | printf("%s", s); |
---|
163 | omFree(s); |
---|
164 | } |
---|
165 | |
---|
166 | bool Fan::addMaxCone(const Cone* pCone) |
---|
167 | { |
---|
168 | /* the given cone must be based on this fan, i.e. refer to |
---|
169 | the same sets of maximal rays and facet normals, etc. */ |
---|
170 | if (pCone->getFan() != this) |
---|
171 | return false; |
---|
172 | |
---|
173 | /* append the adjacency matrix by one more row and column */ |
---|
174 | insertCrosses(1); |
---|
175 | /* in order to use lInsert0 from lists.h/.cc, we pack the cone |
---|
176 | in a list with one element */ |
---|
177 | lists L = (lists)omAllocBin(slists_bin); |
---|
178 | L->Init(1); L->m[0].rtyp = CONE_CMD; L->m[0].data = (void*)pCone; |
---|
179 | /* lInsert0 will create and insert A COPY of its 2nd argument */ |
---|
180 | m_maxCones = lInsert0(m_maxCones, &(L->m[0]), m_maxCones->nr + 1); |
---|
181 | omFreeSize((ADDRESS)L->m, sizeof(sleftv)); |
---|
182 | omFreeBin((ADDRESS)L, slists_bin); |
---|
183 | return true; |
---|
184 | } |
---|
185 | |
---|
186 | int Fan::addMaxCones(const lists L) |
---|
187 | { |
---|
188 | /* the given cones must be based on this fan, i.e. refer to |
---|
189 | the same sets of maximal rays and facet normals, etc. */ |
---|
190 | for (int i = 0; i < L->nr + 1; i++) |
---|
191 | if (((Cone*)(L->m[i].data))->getFan() != this) return i + 1; |
---|
192 | |
---|
193 | /* append the adjacency matrix by sufficiently many new rows and col.'s */ |
---|
194 | insertCrosses(L->nr + 1); |
---|
195 | lists M = lCopy(m_maxCones); |
---|
196 | if (m_maxCones->nr + 1 > 0) |
---|
197 | omFreeSize((ADDRESS)m_maxCones->m, |
---|
198 | (m_maxCones->nr + 1)*sizeof(sleftv)); |
---|
199 | omFreeBin((ADDRESS)m_maxCones, slists_bin); |
---|
200 | m_maxCones = (lists)omAllocBin(slists_bin); |
---|
201 | m_maxCones->Init(M->nr + L->nr + 2); |
---|
202 | for (int i = 0; i < M->nr + 1; i++) |
---|
203 | { |
---|
204 | m_maxCones->m[i].rtyp = CONE_CMD; |
---|
205 | m_maxCones->m[i].data = new Cone(*(Cone*)(M->m[i].data)); |
---|
206 | } |
---|
207 | for (int i = 0; i < L->nr + 1; i++) |
---|
208 | { |
---|
209 | m_maxCones->m[M->nr + 1 + i].rtyp = CONE_CMD; |
---|
210 | m_maxCones->m[M->nr + 1 + i].data = new Cone(*(Cone*)(L->m[i].data)); |
---|
211 | } |
---|
212 | if (M->nr + 1 > 0) |
---|
213 | omFreeSize((ADDRESS)M->m, (M->nr + 1)*sizeof(sleftv)); |
---|
214 | omFreeBin((ADDRESS)M, slists_bin); |
---|
215 | return 0; |
---|
216 | } |
---|
217 | |
---|
218 | bool Fan::contains(const intvec* iv, const int i) |
---|
219 | { |
---|
220 | for (int j = 0; j < iv->length(); j++) |
---|
221 | { |
---|
222 | if ((*iv)[j] == i) return true; |
---|
223 | } |
---|
224 | return false; |
---|
225 | } |
---|
226 | |
---|
227 | bool Fan::deleteMaxCone(const int index) |
---|
228 | { |
---|
229 | /* we want to call deleteMaxCones, thus pack the index in an intvec */ |
---|
230 | intvec* iv = new intvec[1]; (*iv)[0] = index; |
---|
231 | bool result = (deleteMaxCones(iv) == 0); |
---|
232 | delete iv; |
---|
233 | return result; |
---|
234 | } |
---|
235 | |
---|
236 | int Fan::deleteMaxCones(const intvec* indices) |
---|
237 | { |
---|
238 | if (indices == NULL) |
---|
239 | { |
---|
240 | /* delete all maximal cones from this fan */ |
---|
241 | if (m_maxCones->nr + 1 > 0) |
---|
242 | omFreeSize((ADDRESS)m_maxCones->m, |
---|
243 | (m_maxCones->nr + 1)*sizeof(sleftv)); |
---|
244 | omFreeBin((ADDRESS)m_maxCones, slists_bin); |
---|
245 | m_maxCones = (lists)omAllocBin(slists_bin); |
---|
246 | m_maxCones->Init(0); |
---|
247 | /* delete adjacency information */ |
---|
248 | if (m_adjacency->nr + 1 > 0) |
---|
249 | omFreeSize((ADDRESS)m_adjacency->m, |
---|
250 | (m_adjacency->nr + 1)*sizeof(sleftv)); |
---|
251 | omFreeBin((ADDRESS)m_adjacency, slists_bin); |
---|
252 | m_adjacency = (lists)omAllocBin(slists_bin); |
---|
253 | m_adjacency->Init(0); |
---|
254 | return 0; |
---|
255 | } |
---|
256 | else |
---|
257 | { |
---|
258 | int index; |
---|
259 | for (int i = 0; i < indices->length(); i++) |
---|
260 | { |
---|
261 | index = (*indices)[i]; |
---|
262 | if ((index < 1) || (index > m_maxCones->nr + 1)) |
---|
263 | return i + 1; |
---|
264 | } |
---|
265 | /* clean up the adjacency matrix */ |
---|
266 | deleteCrosses(indices); |
---|
267 | lists M = lCopy(m_maxCones); |
---|
268 | if (m_maxCones->nr + 1 > 0) |
---|
269 | omFreeSize((ADDRESS)m_maxCones->m, |
---|
270 | (m_maxCones->nr + 1)*sizeof(sleftv)); |
---|
271 | omFreeBin((ADDRESS)m_maxCones, slists_bin); |
---|
272 | m_maxCones = (lists)omAllocBin(slists_bin); |
---|
273 | m_maxCones->Init(M->nr + 1 - indices->length()); int j = 0; |
---|
274 | for (int i = 0; i < M->nr + 1; i++) |
---|
275 | { |
---|
276 | if (!contains(indices, i + 1)) |
---|
277 | { |
---|
278 | m_maxCones->m[j].rtyp = CONE_CMD; |
---|
279 | m_maxCones->m[j].data = new Cone(*(Cone*)(M->m[i].data)); |
---|
280 | j++; |
---|
281 | } |
---|
282 | } |
---|
283 | if (M->nr + 1 > 0) |
---|
284 | omFreeSize((ADDRESS)M->m, (M->nr + 1)*sizeof(sleftv)); |
---|
285 | omFreeBin((ADDRESS)M, slists_bin); |
---|
286 | return 0; |
---|
287 | } |
---|
288 | } |
---|
289 | |
---|
290 | void Fan::insertCrosses(int s) |
---|
291 | { |
---|
292 | int k = m_maxCones->nr + 1; |
---|
293 | lists L = (lists)omAllocBin(slists_bin); |
---|
294 | L->Init((k + s) * (k + s)); int j = 0; |
---|
295 | for (int i = 0; i < k * k; i++) |
---|
296 | { |
---|
297 | L->m[j ].rtyp = CONE_CMD; |
---|
298 | L->m[j++].data = new Cone(*(Cone*)(m_adjacency->m[i].data)); |
---|
299 | if ((i % k) == k - 1) |
---|
300 | { |
---|
301 | for (int q = 0; q < s; q++) |
---|
302 | { |
---|
303 | L->m[j ].rtyp = CONE_CMD; |
---|
304 | L->m[j++].data = Cone::noAdjacency(); |
---|
305 | } |
---|
306 | } |
---|
307 | } |
---|
308 | for (int q = k * (k + s); q < (k + s) * (k + s); q++) |
---|
309 | { |
---|
310 | L->m[j ].rtyp = CONE_CMD; |
---|
311 | L->m[j++].data = Cone::noAdjacency(); |
---|
312 | } |
---|
313 | if (m_adjacency->nr + 1 > 0) |
---|
314 | omFreeSize((ADDRESS)m_adjacency->m, |
---|
315 | (m_adjacency->nr + 1)*sizeof(sleftv)); |
---|
316 | omFreeBin((ADDRESS)m_adjacency, slists_bin); |
---|
317 | m_adjacency = lCopy(L); |
---|
318 | if (L->nr + 1 > 0) |
---|
319 | omFreeSize((ADDRESS)L->m, |
---|
320 | (L->nr + 1)*sizeof(sleftv)); |
---|
321 | omFreeBin((ADDRESS)L, slists_bin); |
---|
322 | } |
---|
323 | |
---|
324 | void Fan::deleteCrosses(const intvec* indices) |
---|
325 | { |
---|
326 | int k = m_maxCones->nr + 1; |
---|
327 | int d = indices->length(); |
---|
328 | lists L = (lists)omAllocBin(slists_bin); |
---|
329 | L->Init((k - d) * (k - d)); |
---|
330 | int j = 0; |
---|
331 | for (int i = 0; i < k * k; i++) |
---|
332 | { |
---|
333 | while (contains(indices, (i / k) + 1) || contains(indices, (i % k) + 1)) |
---|
334 | i++; /* skip this entry: it's either on a row or on a column |
---|
335 | (or on both) with forbidden index */ |
---|
336 | if (i < k * k) |
---|
337 | { |
---|
338 | L->m[j ].rtyp = CONE_CMD; |
---|
339 | L->m[j++].data = new Cone(*(Cone*)(m_adjacency->m[i].data)); |
---|
340 | } |
---|
341 | } |
---|
342 | if (m_adjacency->nr + 1 > 0) |
---|
343 | omFreeSize((ADDRESS)m_adjacency->m, |
---|
344 | (m_adjacency->nr + 1)*sizeof(sleftv)); |
---|
345 | omFreeBin((ADDRESS)m_adjacency, slists_bin); |
---|
346 | m_adjacency = lCopy(L); |
---|
347 | if (L->nr + 1 > 0) |
---|
348 | omFreeSize((ADDRESS)L->m, |
---|
349 | (L->nr + 1)*sizeof(sleftv)); |
---|
350 | omFreeBin((ADDRESS)L, slists_bin); |
---|
351 | } |
---|
352 | |
---|
353 | Cone* Fan::getAdjacencyFacet(int maxCone1, int maxCone2) const |
---|
354 | { |
---|
355 | int k = getMaxCones()->nr + 1; |
---|
356 | Cone* c = (Cone*)m_adjacency->m[k * (maxCone1 - 1) + maxCone2 - 1].data; |
---|
357 | return new Cone(*c); |
---|
358 | } |
---|
359 | |
---|
360 | intvec* Fan::getAdjacency() const |
---|
361 | { |
---|
362 | int k = getMaxCones()->nr + 1; |
---|
363 | intvec* m = new intvec(k, k, 2); |
---|
364 | Cone* facet; |
---|
365 | for (int r = 1; r <= k; r++) |
---|
366 | for (int c = 1; c <= k; c++) |
---|
367 | { |
---|
368 | facet = (Cone*)m_adjacency->m[k * (r - 1) + c - 1].data; |
---|
369 | if (facet->isNoAdj()) IMATELEM(*m, r, c) = 0; |
---|
370 | else if (facet->isNoFacet()) IMATELEM(*m, r, c) = -1; |
---|
371 | else IMATELEM(*m, r, c) = 1; |
---|
372 | } |
---|
373 | return m; |
---|
374 | } |
---|
375 | |
---|
376 | void Fan::addAdjacency(const int i, const int j, const Cone* c) |
---|
377 | { |
---|
378 | int k = getMaxCones()->nr + 1; |
---|
379 | Cone* d = (Cone*)m_adjacency->m[k * (i - 1) + j - 1].data; |
---|
380 | /* delete old cone */ |
---|
381 | delete d; |
---|
382 | d = (Cone*)m_adjacency->m[k * (j - 1) + i - 1].data; |
---|
383 | /* delete old cone */ |
---|
384 | delete d; |
---|
385 | if (c == NULL) |
---|
386 | { /* symmetric fill */ |
---|
387 | m_adjacency->m[k * (i - 1) + j - 1].data = Cone::noFacet(); |
---|
388 | m_adjacency->m[k * (j - 1) + i - 1].data = Cone::noFacet(); |
---|
389 | } |
---|
390 | else |
---|
391 | { /* symmetric fill */ |
---|
392 | m_adjacency->m[k * (i - 1) + j - 1].data = new Cone(*c); |
---|
393 | m_adjacency->m[k * (j - 1) + i - 1].data = new Cone(*c); |
---|
394 | } |
---|
395 | } |
---|
396 | |
---|
397 | void Fan::addAdjacencies(const intvec* iv, const intvec* jv, const lists L) |
---|
398 | { |
---|
399 | int k = getMaxCones()->nr + 1; |
---|
400 | Cone* d; |
---|
401 | for (int q = 0; q < iv->length(); q++) |
---|
402 | { |
---|
403 | d = (Cone*)m_adjacency->m[k * ((*iv)[q] - 1) + (*jv)[q] - 1].data; |
---|
404 | /* delete old cone */ |
---|
405 | delete d; |
---|
406 | d = (Cone*)m_adjacency->m[k * ((*jv)[q] - 1) + (*iv)[q] - 1].data; |
---|
407 | /* delete old cone */ |
---|
408 | delete d; |
---|
409 | if (L->m[q].rtyp == INT_CMD) /* then (int)(long)L->m[q].data == 0 */ |
---|
410 | { /* symmetric fill */ |
---|
411 | m_adjacency->m[k * ((*iv)[q] - 1) + (*jv)[q] - 1].data = Cone::noFacet(); |
---|
412 | m_adjacency->m[k * ((*jv)[q] - 1) + (*iv)[q] - 1].data = Cone::noFacet(); |
---|
413 | } |
---|
414 | else |
---|
415 | { /* symmetric fill */ |
---|
416 | m_adjacency->m[k * ((*iv)[q] - 1) + (*jv)[q] - 1].data |
---|
417 | = new Cone(*(Cone*)L->m[q].data); |
---|
418 | m_adjacency->m[k * ((*jv)[q] - 1) + (*iv)[q] - 1].data |
---|
419 | = new Cone(*(Cone*)L->m[q].data); |
---|
420 | } |
---|
421 | } |
---|
422 | } |
---|
423 | |
---|
424 | #endif |
---|
425 | /* HAVE_FANS */ |
---|