1 | /* |
---|
2 | Compute the Groebner fan of an ideal |
---|
3 | $Author: monerjan $ |
---|
4 | $Date: 2009/11/03 06:57:32 $ |
---|
5 | $Header: /usr/local/Singular/cvsroot/kernel/gfan.cc,v 1.103 2009/11/03 06:57:32 monerjan Exp $ |
---|
6 | $Id$ |
---|
7 | */ |
---|
8 | |
---|
9 | #include "mod2.h" |
---|
10 | |
---|
11 | #ifdef HAVE_GFAN |
---|
12 | |
---|
13 | #include "kstd1.h" |
---|
14 | #include "kutil.h" //ksCreateSpoly |
---|
15 | // #include "intvec.h" |
---|
16 | #include "polys.h" |
---|
17 | #include "ideals.h" |
---|
18 | #include "kmatrix.h" |
---|
19 | #include "GMPrat.h" |
---|
20 | //#include "fast_maps.h" //Mapping of ideals |
---|
21 | // #include "maps.h" |
---|
22 | // #include "ring.h" //apparently not needed |
---|
23 | // #include "structs.h" |
---|
24 | #include "../Singular/lists.h" |
---|
25 | #include "prCopy.h" //Needed for idrCopyR |
---|
26 | #include <iostream> |
---|
27 | // #include <bitset> |
---|
28 | #include <fstream> //read-write cones to files |
---|
29 | // #include <gmp.h> |
---|
30 | #include <string> |
---|
31 | #include <sstream> |
---|
32 | #include <time.h> |
---|
33 | #include <sys/time.h> |
---|
34 | //#include <gmpxx.h> |
---|
35 | |
---|
36 | /*DO NOT REMOVE THIS*/ |
---|
37 | #ifndef GMPRATIONAL |
---|
38 | #define GMPRATIONAL |
---|
39 | #endif |
---|
40 | |
---|
41 | //Hacks for different working places |
---|
42 | #define p800 |
---|
43 | |
---|
44 | #ifdef UNI |
---|
45 | #include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack |
---|
46 | #include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h" |
---|
47 | #endif |
---|
48 | |
---|
49 | #ifdef HOME |
---|
50 | #include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h" |
---|
51 | #include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h" |
---|
52 | #endif |
---|
53 | |
---|
54 | #ifdef ITWM |
---|
55 | #include "/u/slg/monerjan/cddlib/include/setoper.h" |
---|
56 | #include "/u/slg/monerjan/cddlib/include/cdd.h" |
---|
57 | #include "/u/slg/monerjan/cddlib/include/cddmp.h" |
---|
58 | #endif |
---|
59 | |
---|
60 | #ifdef p800 |
---|
61 | #include "../../cddlib/include/setoper.h" |
---|
62 | #include "../../cddlib/include/cdd.h" |
---|
63 | #include "../../cddlib/include/cddmp.h" |
---|
64 | #endif |
---|
65 | |
---|
66 | #ifndef gfan_DEBUG |
---|
67 | // #define gfan_DEBUG |
---|
68 | #ifndef gfan_DEBUGLEVEL |
---|
69 | #define gfan_DEBUGLEVEL 1 |
---|
70 | #endif |
---|
71 | #endif |
---|
72 | #define gfanp |
---|
73 | #include <gfan.h> |
---|
74 | using namespace std; |
---|
75 | |
---|
76 | /** |
---|
77 | *\brief Class facet |
---|
78 | * Implements the facet structure as a linked list |
---|
79 | * |
---|
80 | */ |
---|
81 | |
---|
82 | /** \brief The default constructor for facets |
---|
83 | */ |
---|
84 | facet::facet() |
---|
85 | { |
---|
86 | // Pointer to next facet. */ |
---|
87 | /* Defaults to NULL. This way there is no need to check explicitly */ |
---|
88 | this->fNormal=NULL; |
---|
89 | this->interiorPoint=NULL; |
---|
90 | this->UCN=0; |
---|
91 | this->codim2Ptr=NULL; |
---|
92 | this->codim=1; //default for (codim-1)-facets |
---|
93 | this->numCodim2Facets=0; |
---|
94 | this->flipGB=NULL; |
---|
95 | this->isIncoming=FALSE; |
---|
96 | this->next=NULL; |
---|
97 | this->prev=NULL; |
---|
98 | this->flipRing=NULL; //the ring on the other side |
---|
99 | this->isFlippable=FALSE; |
---|
100 | } |
---|
101 | |
---|
102 | /** \brief Constructor for facets of codim >= 2 |
---|
103 | * Note that as of now the code of the constructors is only for facets and codim2-faces. One |
---|
104 | * could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar |
---|
105 | */ |
---|
106 | facet::facet(int const &n) |
---|
107 | { |
---|
108 | this->fNormal=NULL; |
---|
109 | this->interiorPoint=NULL; |
---|
110 | this->UCN=0; |
---|
111 | this->codim2Ptr=NULL; |
---|
112 | if(n>1) |
---|
113 | { |
---|
114 | this->codim=n; |
---|
115 | }//NOTE Handle exception here! |
---|
116 | this->numCodim2Facets=0; |
---|
117 | this->flipGB=NULL; |
---|
118 | this->isIncoming=FALSE; |
---|
119 | this->next=NULL; |
---|
120 | this->prev=NULL; |
---|
121 | this->flipRing=NULL; |
---|
122 | this->isFlippable=FALSE; |
---|
123 | } |
---|
124 | |
---|
125 | /** \brief The copy constructor |
---|
126 | * By default only copies the fNormal, f2Normals and UCN |
---|
127 | */ |
---|
128 | facet::facet(const facet& f) |
---|
129 | { |
---|
130 | this->fNormal=ivCopy(f.fNormal); |
---|
131 | this->UCN=f.UCN; |
---|
132 | this->isFlippable=f.isFlippable; |
---|
133 | facet* f2Copy; |
---|
134 | f2Copy=f.codim2Ptr; |
---|
135 | facet* f2Act; |
---|
136 | f2Act=this->codim2Ptr; |
---|
137 | while(f2Copy!=NULL) |
---|
138 | { |
---|
139 | if(f2Act==NULL) |
---|
140 | { |
---|
141 | f2Act=new facet(2); |
---|
142 | this->codim2Ptr=f2Act; |
---|
143 | } |
---|
144 | else |
---|
145 | { |
---|
146 | facet* marker; |
---|
147 | marker = f2Act; |
---|
148 | f2Act->next = new facet(2); |
---|
149 | f2Act = f2Act->next; |
---|
150 | f2Act->prev = marker; |
---|
151 | } |
---|
152 | intvec *f2Normal; |
---|
153 | f2Normal = f2Copy->getFacetNormal(); |
---|
154 | // f2Act->setFacetNormal(f2Copy->getFacetNormal()); |
---|
155 | f2Act->setFacetNormal(f2Normal); |
---|
156 | delete f2Normal; |
---|
157 | f2Act->setUCN(f2Copy->getUCN()); |
---|
158 | f2Copy = f2Copy->next; |
---|
159 | } |
---|
160 | } |
---|
161 | |
---|
162 | /** The default destructor */ |
---|
163 | facet::~facet() |
---|
164 | { |
---|
165 | if(this->fNormal!=NULL) |
---|
166 | delete this->fNormal; |
---|
167 | if(this->interiorPoint!=NULL) |
---|
168 | delete this->interiorPoint; |
---|
169 | /* Cleanup the codim2-structure */ |
---|
170 | if(this->codim==2) |
---|
171 | { |
---|
172 | facet *codim2Ptr; |
---|
173 | codim2Ptr = this->codim2Ptr; |
---|
174 | while(codim2Ptr!=NULL) |
---|
175 | { |
---|
176 | if(codim2Ptr->fNormal!=NULL) |
---|
177 | { |
---|
178 | delete codim2Ptr->fNormal; |
---|
179 | codim2Ptr = codim2Ptr->next; |
---|
180 | } |
---|
181 | } |
---|
182 | } |
---|
183 | if(this->codim2Ptr!=NULL) |
---|
184 | delete this->codim2Ptr; |
---|
185 | if(this->flipGB!=NULL) |
---|
186 | idDelete((ideal *)&this->flipGB); |
---|
187 | if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb) |
---|
188 | // rDelete(this->flipRing); |
---|
189 | // this->flipRing=NULL; |
---|
190 | this->prev=NULL; |
---|
191 | this->next=NULL; |
---|
192 | //this=NULL; |
---|
193 | } |
---|
194 | |
---|
195 | |
---|
196 | /** \brief Comparison of facets*/ |
---|
197 | inline bool gcone::areEqual(facet *f, facet *s) |
---|
198 | { |
---|
199 | bool res = TRUE; |
---|
200 | int notParallelCtr=0; |
---|
201 | int ctr=0; |
---|
202 | intvec* fNormal; //No new since ivCopy and therefore getFacetNormal return a new |
---|
203 | intvec* sNormal; |
---|
204 | fNormal = f->getFacetNormal(); |
---|
205 | sNormal = s->getFacetNormal(); |
---|
206 | if(!isParallel(fNormal,sNormal)) |
---|
207 | notParallelCtr++; |
---|
208 | else//parallelity, so we check the codim2-facets |
---|
209 | { |
---|
210 | facet* f2Act; |
---|
211 | facet* s2Act; |
---|
212 | f2Act = f->codim2Ptr; |
---|
213 | ctr=0; |
---|
214 | while(f2Act!=NULL) |
---|
215 | { |
---|
216 | intvec* f2Normal; |
---|
217 | f2Normal = f2Act->getFacetNormal(); |
---|
218 | s2Act = s->codim2Ptr; |
---|
219 | while(s2Act!=NULL) |
---|
220 | { |
---|
221 | intvec* s2Normal; |
---|
222 | s2Normal = s2Act->getFacetNormal(); |
---|
223 | bool foo=areEqual(f2Normal,s2Normal); |
---|
224 | if(foo) |
---|
225 | ctr++; |
---|
226 | s2Act = s2Act->next; |
---|
227 | delete s2Normal; |
---|
228 | } |
---|
229 | delete f2Normal; |
---|
230 | f2Act = f2Act->next; |
---|
231 | } |
---|
232 | } |
---|
233 | delete fNormal; |
---|
234 | delete sNormal; |
---|
235 | if(ctr==f->numCodim2Facets) |
---|
236 | res=TRUE; |
---|
237 | else |
---|
238 | res=FALSE; |
---|
239 | return res; |
---|
240 | } |
---|
241 | |
---|
242 | /** Stores the facet normal \param intvec*/ |
---|
243 | inline void facet::setFacetNormal(intvec *iv) |
---|
244 | { |
---|
245 | this->fNormal = ivCopy(iv); |
---|
246 | } |
---|
247 | |
---|
248 | /** Hopefully returns the facet normal */ |
---|
249 | inline intvec *facet::getFacetNormal() |
---|
250 | { |
---|
251 | return ivCopy(this->fNormal); |
---|
252 | // return this->fNormal; |
---|
253 | } |
---|
254 | |
---|
255 | /** Method to print the facet normal*/ |
---|
256 | inline void facet::printNormal() |
---|
257 | { |
---|
258 | fNormal->show(); |
---|
259 | } |
---|
260 | |
---|
261 | /** Store the flipped GB*/ |
---|
262 | inline void facet::setFlipGB(ideal I) |
---|
263 | { |
---|
264 | this->flipGB=idCopy(I); |
---|
265 | } |
---|
266 | |
---|
267 | /** Return the flipped GB |
---|
268 | Seems not be used |
---|
269 | Anyhow idCopy would make sense here. |
---|
270 | */ |
---|
271 | inline ideal facet::getFlipGB() |
---|
272 | { |
---|
273 | return this->flipGB; |
---|
274 | } |
---|
275 | |
---|
276 | /** Print the flipped GB*/ |
---|
277 | inline void facet::printFlipGB() |
---|
278 | { |
---|
279 | #ifndef NDEBUG |
---|
280 | idShow(this->flipGB); |
---|
281 | #endif |
---|
282 | } |
---|
283 | |
---|
284 | /** Set the UCN */ |
---|
285 | inline void facet::setUCN(int n) |
---|
286 | { |
---|
287 | this->UCN=n; |
---|
288 | } |
---|
289 | |
---|
290 | /** \brief Get the UCN |
---|
291 | * Returns the UCN iff this != NULL, else -1 |
---|
292 | */ |
---|
293 | inline int facet::getUCN() |
---|
294 | { |
---|
295 | if(this!=NULL)// && this!=(facet *)0xfbfbfbfbfbfbfbfb)// || this!=(facet *)0xfbfbfbfb) ) |
---|
296 | // if(this!=NULL && ( this->fNormal!=(intvec *)0xfbfbfbfb || this->fNormal!=(intvec *)0xfbfbfbfbfbfbfbfb) ) |
---|
297 | return this->UCN; |
---|
298 | else |
---|
299 | return -1; |
---|
300 | } |
---|
301 | |
---|
302 | /** Store an interior point of the facet */ |
---|
303 | inline void facet::setInteriorPoint(intvec *iv) |
---|
304 | { |
---|
305 | this->interiorPoint = ivCopy(iv); |
---|
306 | } |
---|
307 | |
---|
308 | inline intvec *facet::getInteriorPoint() |
---|
309 | { |
---|
310 | return this->interiorPoint; |
---|
311 | } |
---|
312 | |
---|
313 | /** \brief Debugging function |
---|
314 | * prints the facet normal an all (codim-2)-facets that belong to it |
---|
315 | */ |
---|
316 | inline void facet::fDebugPrint() |
---|
317 | { |
---|
318 | facet *codim2Act; |
---|
319 | codim2Act = this->codim2Ptr; |
---|
320 | intvec *fNormal; |
---|
321 | fNormal = this->getFacetNormal(); |
---|
322 | cout << "=======================" << endl; |
---|
323 | cout << "Facet normal = ("; |
---|
324 | fNormal->show(1,1); |
---|
325 | cout << ")"<<endl; |
---|
326 | cout << "-----------------------" << endl; |
---|
327 | cout << "Codim2 facets:" << endl; |
---|
328 | while(codim2Act!=NULL) |
---|
329 | { |
---|
330 | intvec *f2Normal; |
---|
331 | f2Normal = codim2Act->getFacetNormal(); |
---|
332 | cout << "("; |
---|
333 | f2Normal->show(1,0); |
---|
334 | cout << ")" << endl; |
---|
335 | codim2Act = codim2Act->next; |
---|
336 | delete f2Normal; |
---|
337 | } |
---|
338 | cout << "=======================" << endl; |
---|
339 | delete fNormal; |
---|
340 | } |
---|
341 | |
---|
342 | //friend class gcone; //Bad style |
---|
343 | |
---|
344 | |
---|
345 | /** |
---|
346 | *\brief Implements the cone structure |
---|
347 | * |
---|
348 | * A cone is represented by a linked list of facet normals |
---|
349 | * @see facet |
---|
350 | */ |
---|
351 | |
---|
352 | |
---|
353 | /** \brief Default constructor. |
---|
354 | * |
---|
355 | * Initialises this->next=NULL and this->facetPtr=NULL |
---|
356 | */ |
---|
357 | gcone::gcone() |
---|
358 | { |
---|
359 | this->next=NULL; |
---|
360 | this->prev=NULL; |
---|
361 | this->facetPtr=NULL; //maybe this->facetPtr = new facet(); |
---|
362 | this->baseRing=currRing; |
---|
363 | this->counter++; |
---|
364 | this->UCN=this->counter; |
---|
365 | this->numFacets=0; |
---|
366 | this->ivIntPt=NULL; |
---|
367 | } |
---|
368 | |
---|
369 | /** \brief Constructor with ring and ideal |
---|
370 | * |
---|
371 | * This constructor takes the root ring and the root ideal as parameters and stores |
---|
372 | * them in the private members gcone::rootRing and gcone::inputIdeal |
---|
373 | * This constructor is only called once in the computation of the Gröbner fan, |
---|
374 | * namely for the very first cone. Therefore pred is set to 1. |
---|
375 | * Might set it to this->UCN though... |
---|
376 | * Since knowledge of the root ring is only needed when using reverse search, |
---|
377 | * this constructor is not needed when using the "second" method |
---|
378 | */ |
---|
379 | gcone::gcone(ring r, ideal I) |
---|
380 | { |
---|
381 | this->next=NULL; |
---|
382 | this->prev=NULL; |
---|
383 | this->facetPtr=NULL; |
---|
384 | // this->rootRing=r; |
---|
385 | this->inputIdeal=I; |
---|
386 | this->baseRing=currRing; |
---|
387 | this->counter++; |
---|
388 | this->UCN=this->counter; |
---|
389 | this->pred=1; |
---|
390 | this->numFacets=0; |
---|
391 | this->ivIntPt=NULL; |
---|
392 | } |
---|
393 | |
---|
394 | /** \brief Copy constructor |
---|
395 | * |
---|
396 | * Copies a cone, sets this->gcBasis to the flipped GB |
---|
397 | * Call this only after a successful call to gcone::flip which sets facet::flipGB |
---|
398 | */ |
---|
399 | gcone::gcone(const gcone& gc, const facet &f) |
---|
400 | { |
---|
401 | this->next=NULL; |
---|
402 | // this->prev=(gcone *)&gc; //comment in to get a tree |
---|
403 | this->prev=NULL; |
---|
404 | this->numVars=gc.numVars; |
---|
405 | this->counter++; |
---|
406 | this->UCN=this->counter; |
---|
407 | this->pred=gc.UCN; |
---|
408 | this->facetPtr=NULL; |
---|
409 | // this->gcBasis=idCopy(f.flipGB); |
---|
410 | // ring saveRing=currRing; |
---|
411 | // ring tmpRing=f.flipRing; |
---|
412 | // rComplete(f.flipRing); |
---|
413 | // rChangeCurrRing(f.flipRing); |
---|
414 | this->gcBasis=idrCopyR(f.flipGB, f.flipRing); |
---|
415 | // this->inputIdeal=idCopy(this->gcBasis); |
---|
416 | this->baseRing=rCopy(f.flipRing); |
---|
417 | // rComplete(this->baseRing); |
---|
418 | // rChangeCurrRing(saveRing); |
---|
419 | this->numFacets=0; |
---|
420 | this->ivIntPt=NULL; |
---|
421 | // this->rootRing=NULL; |
---|
422 | //rComplete(this->baseRing); |
---|
423 | //rChangeCurrRing(this->baseRing); |
---|
424 | } |
---|
425 | |
---|
426 | /** \brief Default destructor |
---|
427 | */ |
---|
428 | gcone::~gcone() |
---|
429 | { |
---|
430 | if(this->gcBasis!=NULL) |
---|
431 | idDelete((ideal *)&this->gcBasis); |
---|
432 | // if(this->inputIdeal!=NULL) |
---|
433 | // idDelete((ideal *)&this->inputIdeal); |
---|
434 | // if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe) |
---|
435 | // rDelete(this->rootRing); |
---|
436 | // if(this->UCN!=1) |
---|
437 | // rDelete(this->baseRing); |
---|
438 | facet *fAct; |
---|
439 | facet *fDel; |
---|
440 | /*Delete the facet structure*/ |
---|
441 | fAct=this->facetPtr; |
---|
442 | fDel=fAct; |
---|
443 | while(fAct!=NULL) |
---|
444 | { |
---|
445 | fDel=fAct; |
---|
446 | fAct=fAct->next; |
---|
447 | delete fDel; |
---|
448 | } |
---|
449 | this->counter--; |
---|
450 | //should be deleted in noRevS |
---|
451 | // dd_FreeMatrix(this->ddFacets); |
---|
452 | //dd_FreeMatrix(this->ddFacets); |
---|
453 | } |
---|
454 | |
---|
455 | inline int gcone::getCounter() |
---|
456 | { |
---|
457 | return this->counter; |
---|
458 | } |
---|
459 | |
---|
460 | /** \brief Set the interior point of a cone */ |
---|
461 | inline void gcone::setIntPoint(intvec *iv) |
---|
462 | { |
---|
463 | this->ivIntPt=ivCopy(iv); |
---|
464 | } |
---|
465 | |
---|
466 | /** \brief Return the interior point */ |
---|
467 | inline intvec *gcone::getIntPoint() |
---|
468 | { |
---|
469 | return ivCopy(this->ivIntPt); |
---|
470 | } |
---|
471 | |
---|
472 | /** \brief Print the interior point */ |
---|
473 | inline void gcone::showIntPoint() |
---|
474 | { |
---|
475 | ivIntPt->show(); |
---|
476 | } |
---|
477 | |
---|
478 | /** \brief Print facets |
---|
479 | * This is mainly for debugging purposes. Usually called from within gdb |
---|
480 | */ |
---|
481 | inline void gcone::showFacets(const short codim) |
---|
482 | { |
---|
483 | facet *f; |
---|
484 | switch(codim) |
---|
485 | { |
---|
486 | case 1: |
---|
487 | f = this->facetPtr; |
---|
488 | break; |
---|
489 | case 2: |
---|
490 | f = this->facetPtr->codim2Ptr; |
---|
491 | break; |
---|
492 | } |
---|
493 | while(f!=NULL) |
---|
494 | { |
---|
495 | intvec *iv; |
---|
496 | iv = f->getFacetNormal(); |
---|
497 | cout << "("; |
---|
498 | iv->show(1,0); |
---|
499 | if(f->isFlippable==FALSE) |
---|
500 | cout << ")* "; |
---|
501 | else |
---|
502 | cout << ") "; |
---|
503 | delete iv; |
---|
504 | f=f->next; |
---|
505 | } |
---|
506 | cout << endl; |
---|
507 | } |
---|
508 | |
---|
509 | /** For debugging purposes only */ |
---|
510 | inline volatile void gcone::showSLA(facet &f) |
---|
511 | { |
---|
512 | facet *fAct; |
---|
513 | fAct = &f; |
---|
514 | facet *codim2Act; |
---|
515 | codim2Act = fAct->codim2Ptr; |
---|
516 | |
---|
517 | cout << endl; |
---|
518 | while(fAct!=NULL) |
---|
519 | { |
---|
520 | intvec *fNormal; |
---|
521 | fNormal=fAct->getFacetNormal(); |
---|
522 | cout << "("; |
---|
523 | fNormal->show(1,0); |
---|
524 | if(fAct->isFlippable==TRUE) |
---|
525 | cout << ") "; |
---|
526 | else |
---|
527 | cout << ")* "; |
---|
528 | delete fNormal; |
---|
529 | codim2Act = fAct->codim2Ptr; |
---|
530 | cout << " Codim2: "; |
---|
531 | while(codim2Act!=NULL) |
---|
532 | { |
---|
533 | intvec *f2Normal; |
---|
534 | f2Normal = codim2Act->getFacetNormal(); |
---|
535 | cout << "("; |
---|
536 | f2Normal->show(1,0); |
---|
537 | cout << ") "; |
---|
538 | delete f2Normal; |
---|
539 | codim2Act = codim2Act->next; |
---|
540 | } |
---|
541 | cout << "UCN = " << fAct->getUCN() << endl; |
---|
542 | fAct = fAct->next; |
---|
543 | } |
---|
544 | } |
---|
545 | |
---|
546 | inline void gcone::idDebugPrint(const ideal &I) |
---|
547 | { |
---|
548 | int numElts=IDELEMS(I); |
---|
549 | cout << "Ideal with " << numElts << " generators" << endl; |
---|
550 | cout << "Leading terms: "; |
---|
551 | for (int ii=0;ii<numElts;ii++) |
---|
552 | { |
---|
553 | pWrite0(pHead(I->m[ii])); |
---|
554 | cout << ","; |
---|
555 | } |
---|
556 | cout << endl; |
---|
557 | } |
---|
558 | |
---|
559 | inline void gcone::invPrint(const ideal &I) |
---|
560 | { |
---|
561 | // int numElts=IDELEMS(I); |
---|
562 | // cout << "inv = "; |
---|
563 | // for(int ii=0;ii<numElts;ii++); |
---|
564 | // { |
---|
565 | // pWrite0(pHead(I->m[ii])); |
---|
566 | // cout << ","; |
---|
567 | // } |
---|
568 | // cout << endl; |
---|
569 | } |
---|
570 | |
---|
571 | inline bool gcone::isMonomial(const ideal &I) |
---|
572 | { |
---|
573 | bool res = TRUE; |
---|
574 | for(int ii=0;ii<IDELEMS(I);ii++) |
---|
575 | { |
---|
576 | if(pLength((poly)I->m[ii])>1) |
---|
577 | { |
---|
578 | res = FALSE; |
---|
579 | break; |
---|
580 | } |
---|
581 | } |
---|
582 | return res; |
---|
583 | } |
---|
584 | |
---|
585 | /** \brief Set gcone::numFacets */ |
---|
586 | inline void gcone::setNumFacets() |
---|
587 | { |
---|
588 | } |
---|
589 | |
---|
590 | /** \brief Get gcone::numFacets */ |
---|
591 | inline int gcone::getNumFacets() |
---|
592 | { |
---|
593 | return this->numFacets; |
---|
594 | } |
---|
595 | |
---|
596 | inline int gcone::getUCN() |
---|
597 | { |
---|
598 | if( this!=NULL)// && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) ) |
---|
599 | return this->UCN; |
---|
600 | else |
---|
601 | return -1; |
---|
602 | } |
---|
603 | |
---|
604 | inline int gcone::getPredUCN() |
---|
605 | { |
---|
606 | return this->pred; |
---|
607 | } |
---|
608 | |
---|
609 | /** \brief Compute the normals of the cone |
---|
610 | * |
---|
611 | * This method computes a representation of the cone in terms of facet normals. It takes an ideal |
---|
612 | * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize. |
---|
613 | * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44. |
---|
614 | * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects |
---|
615 | * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across |
---|
616 | * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal |
---|
617 | * As a result of this procedure the pointer facetPtr points to the first facet of the cone. |
---|
618 | * |
---|
619 | * Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute |
---|
620 | * an interior point of the cone. |
---|
621 | */ |
---|
622 | inline void gcone::getConeNormals(const ideal &I, bool compIntPoint) |
---|
623 | { |
---|
624 | #ifdef gfanp |
---|
625 | timeval start, end; |
---|
626 | gettimeofday(&start, 0); |
---|
627 | #endif |
---|
628 | poly aktpoly; |
---|
629 | int rows; // will contain the dimensions of the ineq matrix - deprecated by |
---|
630 | dd_rowrange ddrows; |
---|
631 | dd_colrange ddcols; |
---|
632 | dd_rowset ddredrows; // # of redundant rows in ddineq |
---|
633 | dd_rowset ddlinset; // the opposite |
---|
634 | dd_rowindex ddnewpos; // all to make dd_Canonicalize happy |
---|
635 | dd_NumberType ddnumb=dd_Integer; //Number type |
---|
636 | dd_ErrorType dderr=dd_NoError; |
---|
637 | //Compute the # inequalities i.e. rows of the matrix |
---|
638 | rows=0; //Initialization |
---|
639 | for (int ii=0;ii<IDELEMS(I);ii++) |
---|
640 | { |
---|
641 | aktpoly=(poly)I->m[ii]; |
---|
642 | rows=rows+pLength(aktpoly)-1; |
---|
643 | } |
---|
644 | |
---|
645 | dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq |
---|
646 | ddrows=rows; |
---|
647 | ddcols=this->numVars; |
---|
648 | dd_MatrixPtr ddineq; //Matrix to store the inequalities |
---|
649 | ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there |
---|
650 | |
---|
651 | // We loop through each g\in GB and compute the resulting inequalities |
---|
652 | for (int i=0; i<IDELEMS(I); i++) |
---|
653 | { |
---|
654 | aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I |
---|
655 | //simpler version of storing expvect diffs |
---|
656 | int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int)); |
---|
657 | // int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int)); |
---|
658 | pGetExpV(aktpoly,leadexpv); |
---|
659 | while(pNext(aktpoly)!=NULL) |
---|
660 | { |
---|
661 | aktpoly=pNext(aktpoly); |
---|
662 | int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int)); |
---|
663 | pGetExpV(aktpoly,tailexpv); |
---|
664 | for(int kk=1;kk<=this->numVars;kk++) |
---|
665 | { |
---|
666 | dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]); |
---|
667 | } |
---|
668 | aktmatrixrow += 1; |
---|
669 | omFree(tailexpv); |
---|
670 | } |
---|
671 | omFree(leadexpv); |
---|
672 | } //for |
---|
673 | #if true |
---|
674 | /*Let's make the preprocessing here. This could already be done in the above for-loop, |
---|
675 | * but for a start it is more convenient here. |
---|
676 | * We check the necessary condition of FJT p.18 |
---|
677 | * Quote: [...] every non-zero spoly should have at least one of its terms in inv(G) |
---|
678 | */ |
---|
679 | // ideal initialForm=idInit(IDELEMS(I),1); |
---|
680 | intvec *gamma=new intvec(this->numVars); |
---|
681 | int falseGammaCounter=0; |
---|
682 | int *redRowsArray=NULL; |
---|
683 | int num_alloc=0; |
---|
684 | int num_elts=0; |
---|
685 | for(int ii=0;ii<ddineq->rowsize;ii++) |
---|
686 | { |
---|
687 | ideal initialForm=idInit(IDELEMS(I),1); |
---|
688 | //read row ii into gamma |
---|
689 | double tmp; |
---|
690 | for(int jj=1;jj<=this->numVars;jj++) |
---|
691 | { |
---|
692 | tmp=mpq_get_d(ddineq->matrix[ii][jj]); |
---|
693 | (*gamma)[jj-1]=(int)tmp; |
---|
694 | } |
---|
695 | computeInv((ideal&)I,initialForm,*gamma); |
---|
696 | //Create leading ideal |
---|
697 | ideal L=idInit(IDELEMS(initialForm),1); |
---|
698 | for(int jj=0;jj<IDELEMS(initialForm);jj++) |
---|
699 | { |
---|
700 | L->m[jj]=pCopy(pHead(initialForm->m[jj])); |
---|
701 | } |
---|
702 | |
---|
703 | LObject *P = new sLObject(); |
---|
704 | memset(P,0,sizeof(LObject)); |
---|
705 | |
---|
706 | for(int jj=0;jj<=IDELEMS(initialForm)-2;jj++) |
---|
707 | { |
---|
708 | bool isMaybeFacet=FALSE; |
---|
709 | P->p1=initialForm->m[jj]; //build spolys of initialForm in_v |
---|
710 | |
---|
711 | for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++) |
---|
712 | { |
---|
713 | P->p2=initialForm->m[kk]; |
---|
714 | ksCreateSpoly(P); |
---|
715 | if(P->p!=NULL) //spoly non zero=? |
---|
716 | { |
---|
717 | poly p=pInit(); |
---|
718 | poly q=pInit(); |
---|
719 | p=pCopy(P->p); |
---|
720 | q=pHead(p); //Monomial q |
---|
721 | isMaybeFacet=FALSE; |
---|
722 | while(p!=NULL) |
---|
723 | { |
---|
724 | q=pHead(p); |
---|
725 | // unsigned long sevSpoly=pGetShortExpVector(q); |
---|
726 | // unsigned long not_sevL; |
---|
727 | for(int ll=0;ll<IDELEMS(L);ll++) |
---|
728 | { |
---|
729 | // not_sevL=~pGetShortExpVector(L->m[ll]);// |
---|
730 | //if(!(sevSpoly & not_sevL) && pLmDivisibleBy(L->m[ll],q) )//i.e. spoly is in L |
---|
731 | if(pLmEqual(L->m[ll],q) || pDivisibleBy(L->m[ll],q)) |
---|
732 | { |
---|
733 | isMaybeFacet=TRUE; |
---|
734 | break;//for |
---|
735 | } |
---|
736 | } |
---|
737 | if(isMaybeFacet==TRUE) |
---|
738 | { |
---|
739 | break;//while(p!=NULL) |
---|
740 | } |
---|
741 | p=pNext(p); |
---|
742 | }//while |
---|
743 | pDelete(&p); |
---|
744 | pDelete(&q); |
---|
745 | if(isMaybeFacet==FALSE) |
---|
746 | { |
---|
747 | dd_set_si(ddineq->matrix[ii][0],1); |
---|
748 | if(num_alloc==0) |
---|
749 | num_alloc += 1; |
---|
750 | else |
---|
751 | num_alloc += 1; |
---|
752 | void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int))); |
---|
753 | if(!tmp) |
---|
754 | { |
---|
755 | WerrorS("Woah dude! Couldn't realloc memory\n"); |
---|
756 | exit(-1); |
---|
757 | } |
---|
758 | redRowsArray = (int*)tmp; |
---|
759 | redRowsArray[num_elts]=ii; |
---|
760 | num_elts++; |
---|
761 | //break;//for(int kk, since we have found one that is not in L |
---|
762 | goto _start; //mea culpa, mea culpa, mea maxima culpa |
---|
763 | } |
---|
764 | }//if(P->p!=NULL) |
---|
765 | }//for k |
---|
766 | }//for jj |
---|
767 | _start:; |
---|
768 | idDelete(&L); |
---|
769 | delete P; |
---|
770 | idDelete(&initialForm); |
---|
771 | //idDelete(L); |
---|
772 | }//for(ii<ddineq-rowsize |
---|
773 | delete gamma; |
---|
774 | int offset=0;//needed for correction of redRowsArray[ii] |
---|
775 | for( int ii=0;ii<num_elts;ii++ ) |
---|
776 | { |
---|
777 | dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration |
---|
778 | offset++; |
---|
779 | } |
---|
780 | free(redRowsArray); |
---|
781 | /*And now for the strictly positive rows |
---|
782 | * Doesn't gain significant speedup |
---|
783 | */ |
---|
784 | /*int *posRowsArray=NULL; |
---|
785 | num_alloc=0; |
---|
786 | num_elts=0; |
---|
787 | for(int ii=0;ii<ddineq->rowsize;ii++) |
---|
788 | { |
---|
789 | intvec *ivPos = new intvec(this->numVars); |
---|
790 | for(int jj=0;jj<this->numVars;jj++) |
---|
791 | (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]); |
---|
792 | bool isStrictlyPos=FALSE; |
---|
793 | int posCtr=0; |
---|
794 | for(int jj=0;jj<this->numVars;jj++) |
---|
795 | { |
---|
796 | intvec *ivCanonical = new intvec(this->numVars); |
---|
797 | jj==0 ? (*ivCanonical)[ivPos->length()-1]=1 : (*ivCanonical)[jj-1]=1; |
---|
798 | if(dotProduct(*ivCanonical,*ivPos)!=0) |
---|
799 | { |
---|
800 | if ((*ivPos)[jj]>=0) |
---|
801 | { |
---|
802 | posCtr++; |
---|
803 | } |
---|
804 | } |
---|
805 | delete ivCanonical; |
---|
806 | } |
---|
807 | if(posCtr==ivPos->length()) |
---|
808 | isStrictlyPos=TRUE; |
---|
809 | if(isStrictlyPos==TRUE) |
---|
810 | { |
---|
811 | if(num_alloc==0) |
---|
812 | num_alloc += 1; |
---|
813 | else |
---|
814 | num_alloc += 1; |
---|
815 | void *tmp = realloc(posRowsArray,(num_alloc*sizeof(int))); |
---|
816 | if(!tmp) |
---|
817 | { |
---|
818 | WerrorS("Woah dude! Couldn't realloc memory\n"); |
---|
819 | exit(-1); |
---|
820 | } |
---|
821 | posRowsArray = (int*)tmp; |
---|
822 | posRowsArray[num_elts]=ii; |
---|
823 | num_elts++; |
---|
824 | } |
---|
825 | delete ivPos; |
---|
826 | } |
---|
827 | offset=0; |
---|
828 | for(int ii=0;ii<num_elts;ii++) |
---|
829 | { |
---|
830 | dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset); |
---|
831 | offset++; |
---|
832 | } |
---|
833 | free(posRowsArray);*/ |
---|
834 | #endif |
---|
835 | |
---|
836 | dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr); |
---|
837 | //time(&canonicalizeTac); |
---|
838 | //cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl; |
---|
839 | ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed |
---|
840 | ddcols = ddineq->colsize; |
---|
841 | |
---|
842 | //ddCreateMatrix(ddrows,ddcols+1); |
---|
843 | ddFacets = dd_CopyMatrix(ddineq); |
---|
844 | |
---|
845 | /*Write the normals into class facet*/ |
---|
846 | facet *fAct; //pointer to active facet |
---|
847 | int numNonFlip=0; |
---|
848 | for (int kk = 0; kk<ddrows; kk++) |
---|
849 | { |
---|
850 | intvec *load = new intvec(this->numVars); //intvec to store a single facet normal that will then be stored via setFacetNormal |
---|
851 | for (int jj = 1; jj <ddcols; jj++) |
---|
852 | { |
---|
853 | double foo; |
---|
854 | foo = mpq_get_d(ddineq->matrix[kk][jj]); |
---|
855 | (*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load |
---|
856 | }//for (int jj = 1; jj <ddcols; jj++) |
---|
857 | |
---|
858 | /*Quick'n'dirty hack for flippability*/ |
---|
859 | bool isFlip=FALSE; |
---|
860 | for (int jj = 0; jj<load->length(); jj++) |
---|
861 | { |
---|
862 | intvec *ivCanonical = new intvec(load->length()); |
---|
863 | (*ivCanonical)[jj]=1; |
---|
864 | // cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl; |
---|
865 | if (dotProduct(*load,*ivCanonical)<0) |
---|
866 | //if (ivMult(load,ivCanonical)<0) |
---|
867 | { |
---|
868 | isFlip=TRUE; |
---|
869 | break; //URGHS |
---|
870 | } |
---|
871 | delete ivCanonical; |
---|
872 | }/*End of check for flippability*/ |
---|
873 | if (isFlip==FALSE) |
---|
874 | { |
---|
875 | this->numFacets++; |
---|
876 | numNonFlip++; |
---|
877 | if(this->numFacets==1) |
---|
878 | { |
---|
879 | facet *fRoot = new facet(); |
---|
880 | this->facetPtr = fRoot; |
---|
881 | fAct = fRoot; |
---|
882 | |
---|
883 | } |
---|
884 | else |
---|
885 | { |
---|
886 | fAct->next = new facet(); |
---|
887 | fAct = fAct->next; |
---|
888 | } |
---|
889 | fAct->isFlippable=FALSE; |
---|
890 | fAct->setFacetNormal(load); |
---|
891 | fAct->setUCN(this->getUCN()); |
---|
892 | #ifdef gfan_DEBUG |
---|
893 | cout << "Marking facet ("; |
---|
894 | load->show(1,0); |
---|
895 | cout << ") as non flippable" << endl; |
---|
896 | #endif |
---|
897 | } |
---|
898 | else |
---|
899 | { /*Now load should be full and we can call setFacetNormal*/ |
---|
900 | this->numFacets++; |
---|
901 | if(this->numFacets==1) |
---|
902 | { |
---|
903 | facet *fRoot = new facet(); |
---|
904 | this->facetPtr = fRoot; |
---|
905 | fAct = fRoot; |
---|
906 | } |
---|
907 | else |
---|
908 | { |
---|
909 | fAct->next = new facet(); |
---|
910 | fAct = fAct->next; |
---|
911 | } |
---|
912 | fAct->isFlippable=TRUE; |
---|
913 | fAct->setFacetNormal(load); |
---|
914 | fAct->setUCN(this->getUCN()); |
---|
915 | }//if (isFlippable==FALSE) |
---|
916 | delete load; |
---|
917 | }//for (int kk = 0; kk<ddrows; kk++) |
---|
918 | |
---|
919 | //In cases like I=<x-1,y-1> there are only non-flippable facets... |
---|
920 | if(numNonFlip==this->numFacets) |
---|
921 | { |
---|
922 | WerrorS ("Only non-flippable facets. Terminating...\n"); |
---|
923 | // exit(-1);//Bit harsh maybe... |
---|
924 | } |
---|
925 | |
---|
926 | /* |
---|
927 | Now we should have a linked list containing the facet normals of those facets that are |
---|
928 | -irredundant |
---|
929 | -flipable |
---|
930 | Adressing is done via *facetPtr |
---|
931 | */ |
---|
932 | if (compIntPoint==TRUE) |
---|
933 | { |
---|
934 | intvec *iv = new intvec(this->numVars); |
---|
935 | dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1); |
---|
936 | int jj=1; |
---|
937 | for (int ii=0;ii<=this->numVars;ii++) |
---|
938 | { |
---|
939 | dd_set_si(posRestr->matrix[ii][jj],1); |
---|
940 | jj++; |
---|
941 | } |
---|
942 | dd_MatrixAppendTo(&ddineq,posRestr); |
---|
943 | interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets |
---|
944 | this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt |
---|
945 | delete iv; |
---|
946 | dd_FreeMatrix(posRestr); |
---|
947 | } |
---|
948 | //Clean up but don't delete the return value! |
---|
949 | dd_FreeMatrix(ddineq); |
---|
950 | set_free(ddredrows); |
---|
951 | set_free(ddlinset); |
---|
952 | free(ddnewpos); |
---|
953 | #ifdef gfanp |
---|
954 | gettimeofday(&end, 0); |
---|
955 | time_getConeNormals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec)); |
---|
956 | #endif |
---|
957 | |
---|
958 | |
---|
959 | }//gcone::getConeNormals(ideal I) |
---|
960 | |
---|
961 | /** \brief Compute the (codim-2)-facets of a given cone |
---|
962 | * This method is used during noRevS |
---|
963 | * Additionally we check whether the codim2-facet normal is strictly positive. Otherwise |
---|
964 | * the facet is marked as non-flippable. |
---|
965 | */ |
---|
966 | inline void gcone::getCodim2Normals(const gcone &gc) |
---|
967 | { |
---|
968 | #ifdef gfanp |
---|
969 | timeval start, end; |
---|
970 | gettimeofday(&start, 0); |
---|
971 | #endif |
---|
972 | //this->facetPtr->codim2Ptr = new facet(2); //instantiate a (codim-2)-facet |
---|
973 | facet *fAct; |
---|
974 | fAct = this->facetPtr; |
---|
975 | facet *codim2Act; |
---|
976 | //codim2Act = this->facetPtr->codim2Ptr; |
---|
977 | dd_MatrixPtr ddineq;//,P,ddakt; |
---|
978 | dd_ErrorType err; |
---|
979 | //ddineq = facets2Matrix(gc); //get a matrix representation of the cone |
---|
980 | ddineq = dd_CopyMatrix(gc.ddFacets); |
---|
981 | /*Now set appropriate linearity*/ |
---|
982 | for (int ii=0; ii<this->numFacets; ii++) |
---|
983 | { |
---|
984 | // dd_rowset LL; |
---|
985 | // set_initialize(&LL,ddineq->rowsize); |
---|
986 | // set_copy(LL,ddineq->linset); |
---|
987 | // set_copy(ddineq->linset,LL); |
---|
988 | dd_rowset impl_linset, redset; |
---|
989 | dd_rowindex newpos; |
---|
990 | dd_MatrixPtr ddakt; |
---|
991 | ddakt = dd_CopyMatrix(ddineq); |
---|
992 | ddakt->representation=dd_Inequality; |
---|
993 | set_addelem(ddakt->linset,ii+1);/*Now set appropriate linearity*/ |
---|
994 | #ifdef gfanp |
---|
995 | timeval t_ddMC_start, t_ddMC_end; |
---|
996 | gettimeofday(&t_ddMC_start,0); |
---|
997 | #endif |
---|
998 | /*dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);*/ |
---|
999 | // set_copy(LL,ddakt->linset); |
---|
1000 | dd_PolyhedraPtr ddpolyh; |
---|
1001 | ddpolyh=dd_DDMatrix2Poly(ddakt, &err); |
---|
1002 | dd_MatrixPtr P; |
---|
1003 | P=dd_CopyGenerators(ddpolyh); |
---|
1004 | dd_FreePolyhedra(ddpolyh); |
---|
1005 | #ifdef gfanp |
---|
1006 | gettimeofday(&t_ddMC_end,0); |
---|
1007 | t_ddMC += (t_ddMC_end.tv_sec - t_ddMC_start.tv_sec + 1e-6*(t_ddMC_end.tv_usec - t_ddMC_start.tv_usec)); |
---|
1008 | #endif |
---|
1009 | /* We loop through each row of P normalize it by making all |
---|
1010 | * entries integer ones and add the resulting vector to the |
---|
1011 | * int matrix facet::codim2Facets */ |
---|
1012 | for (int jj=1;jj<=/*ddakt*/P->rowsize;jj++) |
---|
1013 | { |
---|
1014 | fAct->numCodim2Facets++; |
---|
1015 | if(fAct->numCodim2Facets==1) |
---|
1016 | { |
---|
1017 | fAct->codim2Ptr = new facet(2); |
---|
1018 | codim2Act = fAct->codim2Ptr; |
---|
1019 | } |
---|
1020 | else |
---|
1021 | { |
---|
1022 | codim2Act->next = new facet(2); |
---|
1023 | codim2Act = codim2Act->next; |
---|
1024 | } |
---|
1025 | intvec *n = new intvec(this->numVars); |
---|
1026 | #ifdef gfanp |
---|
1027 | timeval t_mI_start, t_mI_end; |
---|
1028 | gettimeofday(&t_mI_start,0); |
---|
1029 | #endif |
---|
1030 | makeInt(P,jj,*n); |
---|
1031 | /*for(int kk=0;kk<this->numVars;kk++) |
---|
1032 | { |
---|
1033 | int foo; |
---|
1034 | foo = (int)mpq_get_d(ddakt->matrix[ii][kk+1]); |
---|
1035 | (*n)[kk]=foo; |
---|
1036 | }*/ |
---|
1037 | #ifdef gfanp |
---|
1038 | gettimeofday(&t_mI_end,0); |
---|
1039 | t_mI += (t_mI_end.tv_sec - t_mI_start.tv_sec + 1e-6*(t_mI_end.tv_usec - t_mI_start.tv_usec)); |
---|
1040 | #endif |
---|
1041 | codim2Act->setFacetNormal(n); |
---|
1042 | delete n; |
---|
1043 | } |
---|
1044 | /*We check whether the facet spanned by the codim-2 facets |
---|
1045 | * intersects with the positive orthant. Otherwise we define this |
---|
1046 | * facet to be non-flippable. Works since we set the appropriate |
---|
1047 | * linearity for ddakt above. |
---|
1048 | */ |
---|
1049 | intvec *iv_intPoint = new intvec(this->numVars); |
---|
1050 | dd_MatrixPtr shiftMatrix; |
---|
1051 | dd_MatrixPtr intPointMatrix; |
---|
1052 | shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1); |
---|
1053 | for(int kk=0;kk<this->numVars;kk++) |
---|
1054 | { |
---|
1055 | dd_set_si(shiftMatrix->matrix[kk][0],1); |
---|
1056 | dd_set_si(shiftMatrix->matrix[kk][kk+1],1); |
---|
1057 | } |
---|
1058 | intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix); |
---|
1059 | #ifdef gfanp |
---|
1060 | timeval t_iP_start, t_iP_end; |
---|
1061 | gettimeofday(&t_iP_start, 0); |
---|
1062 | #endif |
---|
1063 | interiorPoint(intPointMatrix,*iv_intPoint); |
---|
1064 | // dd_rowset impl_linste,lbasis; |
---|
1065 | // dd_LPSolutionPtr lps=NULL; |
---|
1066 | // dd_ErrorType err; |
---|
1067 | // dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err); |
---|
1068 | #ifdef gfanp |
---|
1069 | gettimeofday(&t_iP_end, 0); |
---|
1070 | t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec)); |
---|
1071 | #endif |
---|
1072 | for(int ll=0;ll<this->numVars;ll++) |
---|
1073 | { |
---|
1074 | if( (*iv_intPoint)[ll] < 0 ) |
---|
1075 | { |
---|
1076 | fAct->isFlippable=FALSE; |
---|
1077 | break; |
---|
1078 | } |
---|
1079 | } |
---|
1080 | /*End of check*/ |
---|
1081 | fAct = fAct->next; |
---|
1082 | dd_FreeMatrix(ddakt); |
---|
1083 | dd_FreeMatrix(shiftMatrix); |
---|
1084 | dd_FreeMatrix(intPointMatrix); |
---|
1085 | delete iv_intPoint; |
---|
1086 | dd_FreeMatrix(P); |
---|
1087 | // set_free(impl_linset); |
---|
1088 | // set_free(redset); |
---|
1089 | // free(newpos); |
---|
1090 | // set_free(LL); |
---|
1091 | }//for |
---|
1092 | dd_FreeMatrix(ddineq); |
---|
1093 | #ifdef gfanp |
---|
1094 | gettimeofday(&end, 0); |
---|
1095 | time_getCodim2Normals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec)); |
---|
1096 | #endif |
---|
1097 | } |
---|
1098 | |
---|
1099 | /** \brief Compute the Groebner Basis on the other side of a shared facet |
---|
1100 | * |
---|
1101 | * Implements algorithm 4.3.2 from Anders' thesis. |
---|
1102 | * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal |
---|
1103 | * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$ |
---|
1104 | * is parallel to \f$ leadexp(g) \f$ |
---|
1105 | * Parallelity is checked using basic linear algebra. See gcone::isParallel. |
---|
1106 | * Other possibilities include computing the rank of the matrix consisting of the vectors in question and |
---|
1107 | * computing an interior point of the facet and taking all terms having the same weight with respect |
---|
1108 | * to this interior point. |
---|
1109 | *\param ideal, facet |
---|
1110 | * Input: a marked,reduced Groebner basis and a facet |
---|
1111 | */ |
---|
1112 | inline void gcone::flip(ideal gb, facet *f) //Compute "the other side" |
---|
1113 | { |
---|
1114 | #ifdef gfanp |
---|
1115 | timeval start, end; |
---|
1116 | gettimeofday(&start, 0); |
---|
1117 | #endif |
---|
1118 | intvec *fNormal;// = new intvec(this->numVars); //facet normal, check for parallelity |
---|
1119 | fNormal = f->getFacetNormal(); //read this->fNormal; |
---|
1120 | |
---|
1121 | // std::cout << "running gcone::flip" << std::endl; |
---|
1122 | // std::cout << "flipping UCN " << this->getUCN() << endl; |
---|
1123 | // cout << "over facet ("; |
---|
1124 | // fNormal->show(1,0); |
---|
1125 | // cout << ") with UCN " << f->getUCN(); |
---|
1126 | // std::cout << std::endl; |
---|
1127 | if(this->getUCN() != f->getUCN()) |
---|
1128 | { |
---|
1129 | WerrorS("Uh oh... Trying to flip over facet with incompatible UCN"); |
---|
1130 | exit(-1); |
---|
1131 | } |
---|
1132 | /*1st step: Compute the initial ideal*/ |
---|
1133 | /*poly initialFormElement[IDELEMS(gb)];*/ //array of #polys in GB to store initial form |
---|
1134 | ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank); |
---|
1135 | |
---|
1136 | computeInv(gb,initialForm,*fNormal); |
---|
1137 | |
---|
1138 | #ifdef gfan_DEBUG |
---|
1139 | /* cout << "Initial ideal is: " << endl; |
---|
1140 | idShow(initialForm); |
---|
1141 | //f->printFlipGB();*/ |
---|
1142 | // cout << "===" << endl; |
---|
1143 | #endif |
---|
1144 | /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/ |
---|
1145 | /*Substep 2.1 |
---|
1146 | compute $G_{-\alpha}(in_v(I)) |
---|
1147 | see journal p. 66 |
---|
1148 | NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the |
---|
1149 | srcRing already has a weighted ordering |
---|
1150 | */ |
---|
1151 | ring srcRing=currRing; |
---|
1152 | ring tmpRing; |
---|
1153 | |
---|
1154 | if( (srcRing->order[0]!=ringorder_a)) |
---|
1155 | { |
---|
1156 | intvec *iv = new intvec(this->numVars); |
---|
1157 | iv = ivNeg(fNormal); |
---|
1158 | // tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal)); |
---|
1159 | tmpRing=rCopyAndAddWeight(srcRing,iv); |
---|
1160 | delete iv; |
---|
1161 | } |
---|
1162 | else |
---|
1163 | { |
---|
1164 | tmpRing=rCopy0(srcRing); |
---|
1165 | int length=fNormal->length(); |
---|
1166 | int *A=(int *)omAlloc0(length*sizeof(int)); |
---|
1167 | for(int jj=0;jj<length;jj++) |
---|
1168 | { |
---|
1169 | A[jj]=-(*fNormal)[jj]; |
---|
1170 | } |
---|
1171 | omFree(tmpRing->wvhdl[0]); |
---|
1172 | tmpRing->wvhdl[0]=(int*)A; |
---|
1173 | tmpRing->block1[0]=length; |
---|
1174 | rComplete(tmpRing); |
---|
1175 | //omFree(A); |
---|
1176 | } |
---|
1177 | delete fNormal; |
---|
1178 | rChangeCurrRing(tmpRing); |
---|
1179 | |
---|
1180 | ideal ina; |
---|
1181 | ina=idrCopyR(initialForm,srcRing); |
---|
1182 | idDelete(&initialForm); |
---|
1183 | ideal H; |
---|
1184 | // H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous |
---|
1185 | H=kStd(ina,NULL,testHomog,NULL); //This is \mathcal(G)_{>_-\alpha}(in_v(I)) |
---|
1186 | idSkipZeroes(H); |
---|
1187 | idDelete(&ina); |
---|
1188 | |
---|
1189 | /*Substep 2.2 |
---|
1190 | do the lifting and mark according to H |
---|
1191 | */ |
---|
1192 | rChangeCurrRing(srcRing); |
---|
1193 | ideal srcRing_H; |
---|
1194 | ideal srcRing_HH; |
---|
1195 | srcRing_H=idrCopyR(H,tmpRing); |
---|
1196 | //H is needed further below, so don't idDelete here |
---|
1197 | srcRing_HH=ffG(srcRing_H,this->gcBasis); |
---|
1198 | idDelete(&srcRing_H); |
---|
1199 | |
---|
1200 | /*Substep 2.2.1 |
---|
1201 | * Mark according to G_-\alpha |
---|
1202 | * Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis |
---|
1203 | * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone |
---|
1204 | * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha} |
---|
1205 | * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we |
---|
1206 | * compute the difference accordingly |
---|
1207 | */ |
---|
1208 | #ifdef gfanp |
---|
1209 | timeval t_markings_start, t_markings_end; |
---|
1210 | gettimeofday(&t_markings_start, 0); |
---|
1211 | #endif |
---|
1212 | bool markingsAreCorrect=FALSE; |
---|
1213 | dd_MatrixPtr intPointMatrix; |
---|
1214 | int iPMatrixRows=0; |
---|
1215 | dd_rowrange aktrow=0; |
---|
1216 | for (int ii=0;ii<IDELEMS(srcRing_HH);ii++) |
---|
1217 | { |
---|
1218 | poly aktpoly=(poly)srcRing_HH->m[ii];//This is a pointer, so don't pDelete |
---|
1219 | iPMatrixRows = iPMatrixRows+pLength(aktpoly); |
---|
1220 | } |
---|
1221 | /* additionally one row for the standard-simplex and another for a row that becomes 0 during |
---|
1222 | * construction of the differences |
---|
1223 | */ |
---|
1224 | intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); |
---|
1225 | intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational |
---|
1226 | |
---|
1227 | for (int ii=0;ii<IDELEMS(srcRing_HH);ii++) |
---|
1228 | { |
---|
1229 | markingsAreCorrect=FALSE; //crucial to initialise here |
---|
1230 | poly aktpoly=srcRing_HH->m[ii]; //Only a pointer, so don't pDelete |
---|
1231 | /*Comparison of leading monomials is done via exponent vectors*/ |
---|
1232 | for (int jj=0;jj<IDELEMS(H);jj++) |
---|
1233 | { |
---|
1234 | int *src_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int)); |
---|
1235 | int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int)); |
---|
1236 | pGetExpV(aktpoly,src_ExpV); |
---|
1237 | rChangeCurrRing(tmpRing); //this ring change is crucial! |
---|
1238 | pGetExpV(pCopy(H->m[ii]),dst_ExpV); |
---|
1239 | rChangeCurrRing(srcRing); |
---|
1240 | bool expVAreEqual=TRUE; |
---|
1241 | for (int kk=1;kk<=this->numVars;kk++) |
---|
1242 | { |
---|
1243 | #ifdef gfan_DEBUG |
---|
1244 | // cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl; |
---|
1245 | #endif |
---|
1246 | if (src_ExpV[kk]!=dst_ExpV[kk]) |
---|
1247 | { |
---|
1248 | expVAreEqual=FALSE; |
---|
1249 | } |
---|
1250 | } |
---|
1251 | //if (*src_ExpV == *dst_ExpV) |
---|
1252 | if (expVAreEqual==TRUE) |
---|
1253 | { |
---|
1254 | markingsAreCorrect=TRUE; //everything is fine |
---|
1255 | #ifdef gfan_DEBUG |
---|
1256 | // cout << "correct markings" << endl; |
---|
1257 | #endif |
---|
1258 | }//if (pHead(aktpoly)==pHead(H->m[jj]) |
---|
1259 | omFree(src_ExpV); |
---|
1260 | omFree(dst_ExpV); |
---|
1261 | }//for (int jj=0;jj<IDELEMS(H);jj++) |
---|
1262 | |
---|
1263 | int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int)); |
---|
1264 | if (markingsAreCorrect==TRUE) |
---|
1265 | { |
---|
1266 | pGetExpV(aktpoly,leadExpV); |
---|
1267 | } |
---|
1268 | else |
---|
1269 | { |
---|
1270 | rChangeCurrRing(tmpRing); |
---|
1271 | pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial |
---|
1272 | rChangeCurrRing(srcRing); |
---|
1273 | } |
---|
1274 | /*compute differences of the expvects*/ |
---|
1275 | while (pNext(aktpoly)!=NULL) |
---|
1276 | { |
---|
1277 | int *v=(int *)omAlloc((this->numVars+1)*sizeof(int)); |
---|
1278 | /*The following if-else-block makes sure the first term (i.e. the wrongly marked term) |
---|
1279 | is not omitted when computing the differences*/ |
---|
1280 | if(markingsAreCorrect==TRUE) |
---|
1281 | { |
---|
1282 | aktpoly=pNext(aktpoly); |
---|
1283 | pGetExpV(aktpoly,v); |
---|
1284 | } |
---|
1285 | else |
---|
1286 | { |
---|
1287 | pGetExpV(pHead(aktpoly),v); |
---|
1288 | markingsAreCorrect=TRUE; |
---|
1289 | } |
---|
1290 | |
---|
1291 | for (int jj=0;jj<this->numVars;jj++) |
---|
1292 | { |
---|
1293 | /*Store into ddMatrix*/ |
---|
1294 | dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]); |
---|
1295 | } |
---|
1296 | aktrow +=1; |
---|
1297 | omFree(v); |
---|
1298 | } |
---|
1299 | // omFree(v); |
---|
1300 | omFree(leadExpV); |
---|
1301 | }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++) |
---|
1302 | #ifdef gfanp |
---|
1303 | gettimeofday(&t_markings_end, 0); |
---|
1304 | t_markings += (t_markings_end.tv_sec - t_markings_start.tv_sec + 1e-6*(t_markings_end.tv_usec - t_markings_start.tv_usec)); |
---|
1305 | #endif |
---|
1306 | /*Now it is safe to idDelete(H)*/ |
---|
1307 | idDelete(&H); |
---|
1308 | /*Now we add the constraint for the standard simplex*/ |
---|
1309 | // #ifdef gfanp |
---|
1310 | // timeval t_dd_start, t_dd_end; |
---|
1311 | // gettimeofday(&t_dd_start, 0); |
---|
1312 | // #endif |
---|
1313 | dd_set_si(intPointMatrix->matrix[aktrow][0],-1); |
---|
1314 | for (int jj=1;jj<=this->numVars;jj++) |
---|
1315 | { |
---|
1316 | dd_set_si(intPointMatrix->matrix[aktrow][jj],1); |
---|
1317 | } |
---|
1318 | //Let's make sure we compute interior points from the positive orthant |
---|
1319 | dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1); |
---|
1320 | |
---|
1321 | int jj=1; |
---|
1322 | for (int ii=0;ii<this->numVars;ii++) |
---|
1323 | { |
---|
1324 | dd_set_si(posRestr->matrix[ii][jj],1); |
---|
1325 | jj++; |
---|
1326 | } |
---|
1327 | dd_MatrixAppendTo(&intPointMatrix,posRestr); |
---|
1328 | dd_FreeMatrix(posRestr); |
---|
1329 | /*Insert preprocessing here. If it is before the dd_CreateMatrix things go pear shaped |
---|
1330 | * Otherwise things will go pear shaped if called after the standard simplex constraints are added |
---|
1331 | */ |
---|
1332 | preprocessInequalities(intPointMatrix); |
---|
1333 | intvec *iv_weight = new intvec(this->numVars); |
---|
1334 | #ifdef gfanp |
---|
1335 | timeval t_dd_start, t_dd_end; |
---|
1336 | gettimeofday(&t_dd_start, 0); |
---|
1337 | #endif |
---|
1338 | dd_ErrorType err; |
---|
1339 | dd_rowset implLin, redrows; |
---|
1340 | dd_rowindex newpos; |
---|
1341 | |
---|
1342 | dd_MatrixCanonicalize(&intPointMatrix,&implLin,&redrows,&newpos,&err); |
---|
1343 | // dd_MatrixCanonicalizeLinearity(&intPointMatrix,&implLin, &newpos,&err); |
---|
1344 | //dd_MatrixRedundancyRemove is our time sink! |
---|
1345 | // dd_MatrixRedundancyRemove(&intPointMatrix,&redrows,&newpos,&err); |
---|
1346 | interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point |
---|
1347 | dd_FreeMatrix(intPointMatrix); |
---|
1348 | /*Crude attempt for interior point */ |
---|
1349 | /*dd_PolyhedraPtr ddpolyh; |
---|
1350 | dd_ErrorType err; |
---|
1351 | dd_rowset impl_linset,redset; |
---|
1352 | dd_rowindex newpos; |
---|
1353 | dd_MatrixCanonicalize(&intPointMatrix, &impl_linset, &redset, &newpos, &err); |
---|
1354 | ddpolyh=dd_DDMatrix2Poly(intPointMatrix, &err); |
---|
1355 | dd_MatrixPtr P; |
---|
1356 | P=dd_CopyGenerators(ddpolyh); |
---|
1357 | dd_FreePolyhedra(ddpolyh); |
---|
1358 | for(int ii=0;ii<P->rowsize;ii++) |
---|
1359 | { |
---|
1360 | intvec *iv_row=new intvec(this->numVars); |
---|
1361 | makeInt(P,ii+1,*iv_row); |
---|
1362 | iv_weight =ivAdd(iv_weight, iv_row); |
---|
1363 | delete iv_row; |
---|
1364 | } |
---|
1365 | dd_FreeMatrix(P); |
---|
1366 | dd_FreeMatrix(intPointMatrix);*/ |
---|
1367 | #ifdef gfanp |
---|
1368 | gettimeofday(&t_dd_end, 0); |
---|
1369 | t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec)); |
---|
1370 | #endif |
---|
1371 | /*Step 3 |
---|
1372 | turn the minimal basis into a reduced one |
---|
1373 | */ |
---|
1374 | // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a |
---|
1375 | // Thus: |
---|
1376 | //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight); |
---|
1377 | ring dstRing=rCopy0(tmpRing); |
---|
1378 | int length=iv_weight->length(); |
---|
1379 | int *A=(int *)omAlloc0(length*sizeof(int)); |
---|
1380 | for(int jj=0;jj<length;jj++) |
---|
1381 | { |
---|
1382 | A[jj]=(*iv_weight)[jj]; |
---|
1383 | } |
---|
1384 | dstRing->wvhdl[0]=(int*)A; |
---|
1385 | rComplete(dstRing); |
---|
1386 | // omFree(A); |
---|
1387 | rChangeCurrRing(dstRing); |
---|
1388 | rDelete(tmpRing); |
---|
1389 | delete iv_weight; |
---|
1390 | |
---|
1391 | ideal dstRing_I; |
---|
1392 | dstRing_I=idrCopyR(srcRing_HH,srcRing); |
---|
1393 | idDelete(&srcRing_HH); //Hmm.... causes trouble - no more |
---|
1394 | //dstRing_I=idrCopyR(inputIdeal,srcRing); |
---|
1395 | BITSET save=test; |
---|
1396 | test|=Sy_bit(OPT_REDSB); |
---|
1397 | test|=Sy_bit(OPT_REDTAIL); |
---|
1398 | #ifdef gfan_DEBUG |
---|
1399 | test|=Sy_bit(6); //OPT_DEBUG |
---|
1400 | #endif |
---|
1401 | ideal tmpI; |
---|
1402 | //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick |
---|
1403 | //tmpI = idrCopyR(this->inputIdeal,this->baseRing); |
---|
1404 | tmpI = dstRing_I; |
---|
1405 | dstRing_I=kStd(tmpI,NULL,testHomog,NULL); |
---|
1406 | idDelete(&tmpI); |
---|
1407 | idNorm(dstRing_I); |
---|
1408 | // kInterRed(dstRing_I); |
---|
1409 | idSkipZeroes(dstRing_I); |
---|
1410 | test=save; |
---|
1411 | /*End of step 3 - reduction*/ |
---|
1412 | |
---|
1413 | f->setFlipGB(dstRing_I);//store the flipped GB |
---|
1414 | idDelete(&dstRing_I); |
---|
1415 | f->flipRing=rCopy(dstRing); //store the ring on the other side |
---|
1416 | //#ifdef gfan_DEBUG |
---|
1417 | // cout << "Flipped GB is UCN " << counter+1 << ":" << endl; |
---|
1418 | // this->idDebugPrint(dstRing_I); |
---|
1419 | // cout << endl; |
---|
1420 | //#endif |
---|
1421 | rChangeCurrRing(srcRing); //return to the ring we started the computation of flipGB in |
---|
1422 | rDelete(dstRing); |
---|
1423 | #ifdef gfanp |
---|
1424 | gettimeofday(&end, 0); |
---|
1425 | time_flip += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec)); |
---|
1426 | #endif |
---|
1427 | }//void flip(ideal gb, facet *f) |
---|
1428 | |
---|
1429 | /** \brief Compute initial ideal |
---|
1430 | * Compute the initial ideal in_v(G) wrt a (possible) facet normal |
---|
1431 | * used in gcone::getFacetNormal in order to preprocess possible facet normals |
---|
1432 | * and in gcone::flip for obvious reasons. |
---|
1433 | */ |
---|
1434 | inline void gcone::computeInv(ideal &gb, ideal &initialForm, intvec &fNormal) |
---|
1435 | { |
---|
1436 | #ifdef gfanp |
---|
1437 | timeval start, end; |
---|
1438 | gettimeofday(&start, 0); |
---|
1439 | #endif |
---|
1440 | for (int ii=0;ii<IDELEMS(gb);ii++) |
---|
1441 | { |
---|
1442 | poly initialFormElement; |
---|
1443 | poly aktpoly = (poly)gb->m[ii];//Ptr, so don't pDelete(aktpoly) |
---|
1444 | int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int)); |
---|
1445 | pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p) |
---|
1446 | // initialFormElement[ii]=pHead(aktpoly); |
---|
1447 | initialFormElement=pHead(aktpoly); |
---|
1448 | while(pNext(aktpoly)!=NULL) /*loop trough terms and check for parallelity*/ |
---|
1449 | { |
---|
1450 | intvec *check = new intvec(this->numVars); |
---|
1451 | aktpoly=pNext(aktpoly); //next term |
---|
1452 | // pSetm(aktpoly); |
---|
1453 | int *v=(int *)omAlloc((this->numVars+1)*sizeof(int)); |
---|
1454 | pGetExpV(aktpoly,v); |
---|
1455 | /* Convert (int)v into (intvec)check */ |
---|
1456 | for (int jj=0;jj<this->numVars;jj++) |
---|
1457 | { |
---|
1458 | (*check)[jj]=v[jj+1]-leadExpV[jj+1]; |
---|
1459 | } |
---|
1460 | if (isParallel(*check,fNormal)) //pass *check when |
---|
1461 | { |
---|
1462 | //Found a parallel vector. Add it |
---|
1463 | // initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly)); |
---|
1464 | initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly)); |
---|
1465 | } |
---|
1466 | omFree(v); |
---|
1467 | delete check; |
---|
1468 | }//while |
---|
1469 | #ifdef gfan_DEBUG |
---|
1470 | // cout << "Initial Form="; |
---|
1471 | // pWrite(initialFormElement[ii]); |
---|
1472 | // cout << "---" << endl; |
---|
1473 | #endif |
---|
1474 | /*Now initialFormElement must be added to (ideal)initialForm */ |
---|
1475 | // initialForm->m[ii]=pCopy(initialFormElement[ii]); |
---|
1476 | // pDelete(&initialFormElement[ii]); |
---|
1477 | initialForm->m[ii]=pCopy(initialFormElement); |
---|
1478 | pDelete(&initialFormElement); |
---|
1479 | omFree(leadExpV); |
---|
1480 | // delete check; |
---|
1481 | }//for |
---|
1482 | #ifdef gfanp |
---|
1483 | gettimeofday(&end, 0); |
---|
1484 | time_computeInv += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec)); |
---|
1485 | #endif |
---|
1486 | } |
---|
1487 | |
---|
1488 | /** \brief Compute the remainder of a polynomial by a given ideal |
---|
1489 | * |
---|
1490 | * Compute \f$ f^{\mathcal{G}} \f$ |
---|
1491 | * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62 |
---|
1492 | * However, since we are only interested in the remainder, there is no need to |
---|
1493 | * compute the factors \f$ a_i \f$ |
---|
1494 | */ |
---|
1495 | //NOTE: Should be replaced by kNF or kNF2 |
---|
1496 | //NOTE: Done |
---|
1497 | //NOTE: removed with r12286 |
---|
1498 | |
---|
1499 | /** \brief Compute \f$ f-f^{\mathcal{G}} \f$ |
---|
1500 | */ |
---|
1501 | //NOTE: use kNF or kNF2 instead of restOfDivision |
---|
1502 | inline ideal gcone::ffG(const ideal &H, const ideal &G) |
---|
1503 | { |
---|
1504 | // cout << "Entering ffG" << endl; |
---|
1505 | int size=IDELEMS(H); |
---|
1506 | ideal res=idInit(size,1); |
---|
1507 | // poly temp1;//=pInit(); |
---|
1508 | // poly temp2;//=pInit(); |
---|
1509 | // poly temp3;//=pInit(); //polys to temporarily store values for pSub |
---|
1510 | for (int ii=0;ii<size;ii++) |
---|
1511 | { |
---|
1512 | poly temp1=pInit(); |
---|
1513 | poly temp2=pInit(); |
---|
1514 | poly temp3=pInit(); |
---|
1515 | // res->m[ii]=restOfDiv(H->m[ii],G); |
---|
1516 | // res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0)); |
---|
1517 | temp1=pCopy(H->m[ii]); |
---|
1518 | // temp2=pCopy(res->m[ii]); |
---|
1519 | //NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring |
---|
1520 | temp2=pCopy(kNF(G, NULL,H->m[ii],0,0)); |
---|
1521 | temp3=pSub(temp1, temp2); |
---|
1522 | res->m[ii]=pCopy(temp3); |
---|
1523 | //res->m[ii]=pSub(temp1,temp2); //buggy |
---|
1524 | //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]); |
---|
1525 | pDelete(&temp1); |
---|
1526 | // pDelete(&temp2); |
---|
1527 | // pDelete(&temp3); //NOTE does not work, so commented out |
---|
1528 | } |
---|
1529 | return res; |
---|
1530 | } |
---|
1531 | |
---|
1532 | /** \brief Preprocessing of inequlities |
---|
1533 | * Do some preprocessing on the matrix of inequalities |
---|
1534 | * 1) Replace several constraints on the pos. orthants by just one for each orthant |
---|
1535 | * 2) Remove duplicates of inequalities |
---|
1536 | * 3) Remove inequalities that arise as sums of other inequalities |
---|
1537 | */ |
---|
1538 | void gcone::preprocessInequalities(dd_MatrixPtr &ddineq) |
---|
1539 | { |
---|
1540 | //Remove strictly positive rows |
---|
1541 | // int *posRowsArray=NULL; |
---|
1542 | // int num_alloc=0; |
---|
1543 | // int num_elts=0; |
---|
1544 | // for(int ii=0;ii<ddineq->rowsize;ii++) |
---|
1545 | // { |
---|
1546 | // intvec *ivPos = new intvec(pVariables); |
---|
1547 | // for(int jj=0;jj<pVariables;jj++) |
---|
1548 | // (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]); |
---|
1549 | // bool isStrictlyPos=FALSE; |
---|
1550 | // int posCtr=0; |
---|
1551 | // for(int jj=0;jj<pVariables;jj++) |
---|
1552 | // { |
---|
1553 | // intvec *ivCanonical = new intvec(pVariables); |
---|
1554 | // jj==0 ? (*ivCanonical)[ivPos->length()-1]=1 : (*ivCanonical)[jj-1]=1; |
---|
1555 | // if(dotProduct(*ivCanonical,*ivPos)!=0) |
---|
1556 | // { |
---|
1557 | // if ((*ivPos)[jj]>=0) |
---|
1558 | // { |
---|
1559 | // posCtr++; |
---|
1560 | // } |
---|
1561 | // } |
---|
1562 | // delete ivCanonical; |
---|
1563 | // } |
---|
1564 | // if(posCtr==ivPos->length()) |
---|
1565 | // isStrictlyPos=TRUE; |
---|
1566 | // if(isStrictlyPos==TRUE) |
---|
1567 | // { |
---|
1568 | // if(num_alloc==0) |
---|
1569 | // num_alloc += 1; |
---|
1570 | // else |
---|
1571 | // num_alloc += 1; |
---|
1572 | // void *tmp = realloc(posRowsArray,(num_alloc*sizeof(int))); |
---|
1573 | // if(!tmp) |
---|
1574 | // { |
---|
1575 | // WerrorS("Woah dude! Couldn't realloc memory\n"); |
---|
1576 | // exit(-1); |
---|
1577 | // } |
---|
1578 | // posRowsArray = (int*)tmp; |
---|
1579 | // posRowsArray[num_elts]=ii; |
---|
1580 | // num_elts++; |
---|
1581 | // } |
---|
1582 | // delete ivPos; |
---|
1583 | // } |
---|
1584 | int offset=0; |
---|
1585 | // for(int ii=0;ii<num_elts;ii++) |
---|
1586 | // { |
---|
1587 | // dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset); |
---|
1588 | // offset++; |
---|
1589 | // } |
---|
1590 | // free(posRowsArray); |
---|
1591 | //Remove zeroes |
---|
1592 | int rowsize=ddineq->rowsize; |
---|
1593 | for(int ii=0;ii<ddineq->rowsize;ii++) |
---|
1594 | { |
---|
1595 | intvec *iv = new intvec(this->numVars); |
---|
1596 | int posCtr=0; |
---|
1597 | for(int jj=0;jj<this->numVars;jj++) |
---|
1598 | { |
---|
1599 | (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]); |
---|
1600 | if((*iv)[ii]>0) |
---|
1601 | posCtr++; |
---|
1602 | } |
---|
1603 | if( (iv->compare(0)==0)) //|| (posCtr==iv->length()) ) |
---|
1604 | { |
---|
1605 | dd_MatrixRowRemove(&ddineq,ii+1); |
---|
1606 | ii--; |
---|
1607 | // rowsize=ddineq->rowsize; |
---|
1608 | } |
---|
1609 | delete iv; |
---|
1610 | } |
---|
1611 | //Remove duplicates of rows |
---|
1612 | // posRowsArray=NULL; |
---|
1613 | // num_alloc=0; |
---|
1614 | // num_elts=0; |
---|
1615 | // offset=0; |
---|
1616 | // int num_newRows = ddineq->rowsize; |
---|
1617 | // for(int ii=0;ii<ddineq->rowsize-1;ii++) |
---|
1618 | // for(int ii=0;ii<num_newRows-1;ii++) |
---|
1619 | // { |
---|
1620 | // intvec *iv = new intvec(this->numVars);//1st vector to check against |
---|
1621 | // for(int jj=0;jj<this->numVars;jj++) |
---|
1622 | // (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]); |
---|
1623 | // for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++) |
---|
1624 | // { |
---|
1625 | // intvec *ivCheck = new intvec(this->numVars);//Checked against iv |
---|
1626 | // for(int kk=0;kk<this->numVars;kk++) |
---|
1627 | // (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]); |
---|
1628 | // if (iv->compare(ivCheck)==0) |
---|
1629 | // { |
---|
1630 | // // cout << "=" << endl; |
---|
1631 | // // num_alloc++; |
---|
1632 | // // void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int))); |
---|
1633 | // // if(!tmp) |
---|
1634 | // // { |
---|
1635 | // // WerrorS("Woah dude! Couldn't realloc memory\n"); |
---|
1636 | // // exit(-1); |
---|
1637 | // // } |
---|
1638 | // // posRowsArray = (int*)tmp; |
---|
1639 | // // posRowsArray[num_elts]=jj; |
---|
1640 | // // num_elts++; |
---|
1641 | // dd_MatrixRowRemove(&ddineq,jj+1); |
---|
1642 | // num_newRows = ddineq->rowsize; |
---|
1643 | // } |
---|
1644 | // delete ivCheck; |
---|
1645 | // } |
---|
1646 | // delete iv; |
---|
1647 | // } |
---|
1648 | // for(int ii=0;ii<num_elts;ii++) |
---|
1649 | // { |
---|
1650 | // dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset); |
---|
1651 | // offset++; |
---|
1652 | // } |
---|
1653 | // free(posRowsArray); |
---|
1654 | //Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/ |
---|
1655 | }//preprocessInequalities |
---|
1656 | |
---|
1657 | /** \brief Compute a Groebner Basis |
---|
1658 | * |
---|
1659 | * Computes the Groebner basis and stores the result in gcone::gcBasis |
---|
1660 | *\param ideal |
---|
1661 | *\return void |
---|
1662 | */ |
---|
1663 | inline void gcone::getGB(const ideal &inputIdeal) |
---|
1664 | { |
---|
1665 | BITSET save=test; |
---|
1666 | test|=Sy_bit(OPT_REDSB); |
---|
1667 | test|=Sy_bit(OPT_REDTAIL); |
---|
1668 | ideal gb; |
---|
1669 | gb=kStd(inputIdeal,NULL,testHomog,NULL); |
---|
1670 | idNorm(gb); |
---|
1671 | idSkipZeroes(gb); |
---|
1672 | this->gcBasis=gb; //write the GB into gcBasis |
---|
1673 | test=save; |
---|
1674 | }//void getGB |
---|
1675 | |
---|
1676 | /** \brief Compute the negative of a given intvec |
---|
1677 | */ |
---|
1678 | inline intvec *gcone::ivNeg(const intvec *iv) |
---|
1679 | { |
---|
1680 | //NOTE: Can't this be done without new? |
---|
1681 | intvec *res;// = new intvec(iv->length()); |
---|
1682 | res=ivCopy(iv); |
---|
1683 | *res *= (int)-1; |
---|
1684 | return res; |
---|
1685 | } |
---|
1686 | |
---|
1687 | |
---|
1688 | /** \brief Compute the dot product of two intvecs |
---|
1689 | * |
---|
1690 | */ |
---|
1691 | inline int gcone::dotProduct(const intvec &iva, const intvec &ivb) |
---|
1692 | { |
---|
1693 | int res=0; |
---|
1694 | for (int i=0;i<this->numVars;i++) |
---|
1695 | { |
---|
1696 | res = res+(iva[i]*ivb[i]); |
---|
1697 | } |
---|
1698 | return res; |
---|
1699 | }//int dotProduct |
---|
1700 | |
---|
1701 | /** \brief Check whether two intvecs are parallel |
---|
1702 | * |
---|
1703 | * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$ |
---|
1704 | */ |
---|
1705 | inline bool gcone::isParallel(const intvec &a, const intvec &b) |
---|
1706 | { |
---|
1707 | int lhs,rhs; |
---|
1708 | bool res; |
---|
1709 | lhs=dotProduct(a,b)*dotProduct(a,b); |
---|
1710 | rhs=dotProduct(a,a)*dotProduct(b,b); |
---|
1711 | //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl; |
---|
1712 | if (lhs==rhs) |
---|
1713 | { |
---|
1714 | res = TRUE; |
---|
1715 | } |
---|
1716 | else |
---|
1717 | { |
---|
1718 | res = FALSE; |
---|
1719 | } |
---|
1720 | return res; |
---|
1721 | }//bool isParallel |
---|
1722 | |
---|
1723 | /** \brief Compute an interior point of a given cone |
---|
1724 | * Result will be written into intvec iv. |
---|
1725 | * Any rational point is automatically converted into an integer. |
---|
1726 | */ |
---|
1727 | inline void gcone::interiorPoint( dd_MatrixPtr &M, intvec &iv) //no const &M here since we want to remove redundant rows |
---|
1728 | { |
---|
1729 | dd_LPPtr lp,lpInt; |
---|
1730 | dd_ErrorType err=dd_NoError; |
---|
1731 | dd_LPSolverType solver=dd_DualSimplex; |
---|
1732 | dd_LPSolutionPtr lpSol=NULL; |
---|
1733 | dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior |
---|
1734 | dd_rowindex ddnewpos; |
---|
1735 | dd_NumberType numb; |
---|
1736 | //M->representation=dd_Inequality; |
---|
1737 | |
---|
1738 | //NOTE: Make this n-dimensional! |
---|
1739 | //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1); |
---|
1740 | |
---|
1741 | /*NOTE: Leave the following line commented out! |
---|
1742 | * Otherwise it will slow down computations a great deal |
---|
1743 | * */ |
---|
1744 | // dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err); |
---|
1745 | //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;} |
---|
1746 | dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1); |
---|
1747 | int jj=1; |
---|
1748 | for (int ii=0;ii<this->numVars;ii++) |
---|
1749 | { |
---|
1750 | dd_set_si(posRestr->matrix[ii][jj],1); |
---|
1751 | jj++; |
---|
1752 | } |
---|
1753 | dd_MatrixAppendTo(&M,posRestr); |
---|
1754 | dd_FreeMatrix(posRestr); |
---|
1755 | // dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err); |
---|
1756 | lp=dd_Matrix2LP(M, &err); |
---|
1757 | if (err!=dd_NoError){WerrorS("Error during dd_Matrix2LP in gcone::interiorPoint");} |
---|
1758 | if (lp==NULL){WerrorS("LP is NULL");} |
---|
1759 | #ifdef gfan_DEBUG |
---|
1760 | // dd_WriteLP(stdout,lp); |
---|
1761 | #endif |
---|
1762 | |
---|
1763 | lpInt=dd_MakeLPforInteriorFinding(lp); |
---|
1764 | if (err!=dd_NoError){WerrorS("Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint");} |
---|
1765 | #ifdef gfan_DEBUG |
---|
1766 | // dd_WriteLP(stdout,lpInt); |
---|
1767 | #endif |
---|
1768 | |
---|
1769 | // dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err); |
---|
1770 | if (err!=dd_NoError) |
---|
1771 | { |
---|
1772 | WerrorS("Error during dd_FindRelativeInterior in gcone::interiorPoint"); |
---|
1773 | dd_WriteErrorMessages(stdout, err); |
---|
1774 | } |
---|
1775 | |
---|
1776 | dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior |
---|
1777 | // if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");} |
---|
1778 | |
---|
1779 | lpSol=dd_CopyLPSolution(lpInt); |
---|
1780 | // if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");} |
---|
1781 | #ifdef gfan_DEBUG |
---|
1782 | cout << "Interior point: "; |
---|
1783 | for (int ii=1; ii<(lpSol->d)-1;ii++) |
---|
1784 | { |
---|
1785 | dd_WriteNumber(stdout,lpSol->sol[ii]); |
---|
1786 | } |
---|
1787 | cout << endl; |
---|
1788 | #endif |
---|
1789 | |
---|
1790 | //NOTE The following strongly resembles parts of makeInt. |
---|
1791 | //Maybe merge sometimes |
---|
1792 | mpz_t kgV; mpz_init(kgV); |
---|
1793 | mpz_set_str(kgV,"1",10); |
---|
1794 | mpz_t den; mpz_init(den); |
---|
1795 | mpz_t tmp; mpz_init(tmp); |
---|
1796 | mpq_get_den(tmp,lpSol->sol[1]); |
---|
1797 | for (int ii=1;ii<(lpSol->d)-1;ii++) |
---|
1798 | { |
---|
1799 | mpq_get_den(den,lpSol->sol[ii+1]); |
---|
1800 | mpz_lcm(kgV,tmp,den); |
---|
1801 | mpz_set(tmp, kgV); |
---|
1802 | } |
---|
1803 | mpq_t qkgV; |
---|
1804 | mpq_init(qkgV); |
---|
1805 | mpq_set_z(qkgV,kgV); |
---|
1806 | for (int ii=1;ii<(lpSol->d)-1;ii++) |
---|
1807 | { |
---|
1808 | mpq_t product; |
---|
1809 | mpq_init(product); |
---|
1810 | mpq_mul(product,qkgV,lpSol->sol[ii]); |
---|
1811 | iv[ii-1]=(int)mpz_get_d(mpq_numref(product)); |
---|
1812 | mpq_clear(product); |
---|
1813 | } |
---|
1814 | #ifdef gfan_DEBUG |
---|
1815 | // iv.show(); |
---|
1816 | // cout << endl; |
---|
1817 | #endif |
---|
1818 | mpq_clear(qkgV); |
---|
1819 | mpz_clear(tmp); |
---|
1820 | mpz_clear(den); |
---|
1821 | mpz_clear(kgV); |
---|
1822 | |
---|
1823 | dd_FreeLPSolution(lpSol); |
---|
1824 | dd_FreeLPData(lpInt); |
---|
1825 | dd_FreeLPData(lp); |
---|
1826 | // set_free(ddlinset); |
---|
1827 | // set_free(ddredrows); |
---|
1828 | |
---|
1829 | }//void interiorPoint(dd_MatrixPtr const &M) |
---|
1830 | |
---|
1831 | inline void gcone::interiorPoint2(const dd_MatrixPtr &M, intvec &iv) |
---|
1832 | { |
---|
1833 | KMatrix<Rational> mMat(M->rowsize+1,M->colsize); |
---|
1834 | for(int ii=0;ii<M->rowsize;ii++) |
---|
1835 | { |
---|
1836 | for(int jj=0;jj<M->colsize-1;jj++) |
---|
1837 | { |
---|
1838 | if(mpq_sgn(M->matrix[ii][jj+1])<-1) |
---|
1839 | { |
---|
1840 | mMat.set(ii,jj,-(Rational)mpq_get_d(M->matrix[ii][jj+1])); |
---|
1841 | } |
---|
1842 | else |
---|
1843 | mMat.set(ii,jj,(Rational)mpq_get_d(M->matrix[ii][jj+1])); |
---|
1844 | |
---|
1845 | // mMat.set(ii,jj,&(M->matrix[ii][jj+1]) ); |
---|
1846 | cout << mpq_get_d(M->matrix[ii][jj+1]) << " "; |
---|
1847 | // int val=(int)mMat.get(ii,jj); |
---|
1848 | // cout << ii << "," << jj << endl;; |
---|
1849 | // mpq_out_str (NULL, 10, (__mpq_struct)mMat.get(ii,jj)); |
---|
1850 | } |
---|
1851 | cout << endl; |
---|
1852 | mMat.set(ii,M->colsize-1,1); |
---|
1853 | } |
---|
1854 | dd_WriteMatrix(stdout,M); |
---|
1855 | // for(int ii=0;ii<M->rowsize;ii++) |
---|
1856 | // { |
---|
1857 | // cout << mMat.get(ii,ii+M->colsize) << " "; |
---|
1858 | // if((ii+M->colsize)%M->colsize==0) |
---|
1859 | // cout << endl; |
---|
1860 | // } |
---|
1861 | |
---|
1862 | Rational* mSol; |
---|
1863 | int rank; |
---|
1864 | int c; |
---|
1865 | // dd_WriteMatrix(stdout,M); |
---|
1866 | rank=mMat.solve(&mSol,&c); |
---|
1867 | // for(int ii=0;ii<c;ii++) |
---|
1868 | // iv[ii]=mSol[ii]; |
---|
1869 | // cout << mSol[ii].get_den() << "/" << mSol[ii].get_num() << endl; |
---|
1870 | int gcd=1; |
---|
1871 | for(int ii=0;ii<c-1;ii++) |
---|
1872 | gcd += intgcd(mSol[ii].get_den_si(),mSol[ii+1].get_den_si()); |
---|
1873 | cout << gcd << endl; |
---|
1874 | for(int ii=0;ii<iv.length();ii++) |
---|
1875 | iv[ii]=(int)((mSol[ii].get_num_si()*gcd)/mSol[ii].get_den_si()); |
---|
1876 | |
---|
1877 | } |
---|
1878 | |
---|
1879 | |
---|
1880 | /** \brief Copy a ring and add a weighted ordering in first place |
---|
1881 | * |
---|
1882 | */ |
---|
1883 | ring gcone::rCopyAndAddWeight(const ring &r, const intvec *ivw) |
---|
1884 | { |
---|
1885 | ring res=rCopy0(r); |
---|
1886 | int jj; |
---|
1887 | |
---|
1888 | omFree(res->order); |
---|
1889 | res->order =(int *)omAlloc0(4*sizeof(int)); |
---|
1890 | omFree(res->block0); |
---|
1891 | res->block0=(int *)omAlloc0(4*sizeof(int)); |
---|
1892 | omFree(res->block1); |
---|
1893 | res->block1=(int *)omAlloc0(4*sizeof(int)); |
---|
1894 | omfree(res->wvhdl); |
---|
1895 | res->wvhdl =(int **)omAlloc0(4*sizeof(int**)); |
---|
1896 | |
---|
1897 | res->order[0]=ringorder_a; |
---|
1898 | res->block0[0]=1; |
---|
1899 | res->block1[0]=res->N;; |
---|
1900 | res->order[1]=ringorder_dp; //basically useless, since that should never be used |
---|
1901 | res->block0[1]=1; |
---|
1902 | res->block1[1]=res->N;; |
---|
1903 | res->order[2]=ringorder_C; |
---|
1904 | |
---|
1905 | int length=ivw->length(); |
---|
1906 | int *A=(int *)omAlloc0(length*sizeof(int)); |
---|
1907 | for (jj=0;jj<length;jj++) |
---|
1908 | { |
---|
1909 | A[jj]=(*ivw)[jj]; |
---|
1910 | } |
---|
1911 | res->wvhdl[0]=(int *)A; |
---|
1912 | res->block1[0]=length; |
---|
1913 | |
---|
1914 | rComplete(res); |
---|
1915 | return res; |
---|
1916 | }//rCopyAndAdd |
---|
1917 | |
---|
1918 | ring rCopyAndChangeWeight(ring const &r, intvec *ivw) |
---|
1919 | { |
---|
1920 | ring res=rCopy0(currRing); |
---|
1921 | rComplete(res); |
---|
1922 | rSetWeightVec(res,(int64*)ivw); |
---|
1923 | //rChangeCurrRing(rnew); |
---|
1924 | return res; |
---|
1925 | } |
---|
1926 | |
---|
1927 | /** \brief Checks whether a given facet is a search facet |
---|
1928 | * Determines whether a given facet of a cone is the search facet of a neighbouring cone |
---|
1929 | * This is done in the following way: |
---|
1930 | * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet |
---|
1931 | * that is first crossed during the generic walk. |
---|
1932 | * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet. |
---|
1933 | * If this is the case, then our facet is indeed a search facet and TRUE is retuned. |
---|
1934 | */ |
---|
1935 | //removed with r12286 |
---|
1936 | |
---|
1937 | /** \brief Check for equality of two intvecs |
---|
1938 | */ |
---|
1939 | inline bool gcone::areEqual(const intvec &a, const intvec &b) |
---|
1940 | { |
---|
1941 | bool res=TRUE; |
---|
1942 | for(int ii=0;ii<this->numVars;ii++) |
---|
1943 | { |
---|
1944 | if(a[ii]!=b[ii]) |
---|
1945 | { |
---|
1946 | res=FALSE; |
---|
1947 | break; |
---|
1948 | } |
---|
1949 | } |
---|
1950 | return res; |
---|
1951 | } |
---|
1952 | |
---|
1953 | /** \brief The reverse search algorithm |
---|
1954 | */ |
---|
1955 | //removed with r12286 |
---|
1956 | |
---|
1957 | /** \brief The new method of Markwig and Jensen |
---|
1958 | * Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals. |
---|
1959 | * In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components. |
---|
1960 | * Keep a list of facets as a linked list containing an intvec and an integer matrix. |
---|
1961 | * Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as |
---|
1962 | * soon as we compute this facet again. Comparison of facets is done by... |
---|
1963 | */ |
---|
1964 | void gcone::noRevS(gcone &gcRoot, bool usingIntPoint) |
---|
1965 | { |
---|
1966 | facet *SearchListRoot = new facet(); //The list containing ALL facets we come across |
---|
1967 | facet *SearchListAct; |
---|
1968 | SearchListAct = NULL; |
---|
1969 | //SearchListAct = SearchListRoot; |
---|
1970 | |
---|
1971 | gcone *gcAct; |
---|
1972 | gcAct = &gcRoot; |
---|
1973 | gcone *gcPtr; //Pointer to end of linked list of cones |
---|
1974 | gcPtr = &gcRoot; |
---|
1975 | gcone *gcNext; //Pointer to next cone to be visited |
---|
1976 | gcNext = &gcRoot; |
---|
1977 | gcone *gcHead; |
---|
1978 | gcHead = &gcRoot; |
---|
1979 | |
---|
1980 | facet *fAct; |
---|
1981 | fAct = gcAct->facetPtr; |
---|
1982 | |
---|
1983 | ring rAct; |
---|
1984 | rAct = currRing; |
---|
1985 | |
---|
1986 | int UCNcounter=gcAct->getUCN(); |
---|
1987 | |
---|
1988 | /* Use pure SLA or writeCone2File */ |
---|
1989 | // enum heuristic { |
---|
1990 | // ram, |
---|
1991 | // hdd |
---|
1992 | // }; |
---|
1993 | // heuristic h; |
---|
1994 | // h=hdd; |
---|
1995 | |
---|
1996 | #ifdef gfan_DEBUG |
---|
1997 | cout << "NoRevs" << endl; |
---|
1998 | cout << "Facets are:" << endl; |
---|
1999 | gcAct->showFacets(); |
---|
2000 | #endif |
---|
2001 | |
---|
2002 | gcAct->getCodim2Normals(*gcAct); |
---|
2003 | |
---|
2004 | //Compute unique representation of codim-2-facets |
---|
2005 | gcAct->normalize(); |
---|
2006 | |
---|
2007 | /* Make a copy of the facet list of first cone |
---|
2008 | Since the operations getCodim2Normals and normalize affect the facets |
---|
2009 | we must not memcpy them before these ops! |
---|
2010 | */ |
---|
2011 | |
---|
2012 | facet *codim2Act; codim2Act = NULL; |
---|
2013 | facet *sl2Root; //sl2Root = new facet(2); |
---|
2014 | facet *sl2Act; //sl2Act = sl2Root; |
---|
2015 | |
---|
2016 | for(int ii=0;ii<this->numFacets;ii++) |
---|
2017 | { |
---|
2018 | //only copy flippable facets! NOTE: We do not need flipRing or any such information. |
---|
2019 | if(fAct->isFlippable==TRUE) |
---|
2020 | { |
---|
2021 | intvec *fNormal; |
---|
2022 | fNormal = fAct->getFacetNormal(); |
---|
2023 | if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable |
---|
2024 | { |
---|
2025 | SearchListAct = SearchListRoot; |
---|
2026 | } |
---|
2027 | else |
---|
2028 | { |
---|
2029 | SearchListAct->next = new facet(); |
---|
2030 | SearchListAct = SearchListAct->next; |
---|
2031 | } |
---|
2032 | SearchListAct->setFacetNormal(fNormal); |
---|
2033 | SearchListAct->setUCN(this->getUCN()); |
---|
2034 | SearchListAct->numCodim2Facets=fAct->numCodim2Facets; |
---|
2035 | SearchListAct->isFlippable=TRUE; |
---|
2036 | //Copy codim2-facets |
---|
2037 | codim2Act=fAct->codim2Ptr; |
---|
2038 | SearchListAct->codim2Ptr = new facet(2); |
---|
2039 | sl2Root = SearchListAct->codim2Ptr; |
---|
2040 | sl2Act = sl2Root; |
---|
2041 | //while(codim2Act!=NULL) |
---|
2042 | for(int jj=0;jj<fAct->numCodim2Facets;jj++) |
---|
2043 | { |
---|
2044 | intvec *f2Normal; |
---|
2045 | f2Normal = codim2Act->getFacetNormal(); |
---|
2046 | if(jj==0) |
---|
2047 | { |
---|
2048 | sl2Act = sl2Root; |
---|
2049 | sl2Act->setFacetNormal(f2Normal); |
---|
2050 | } |
---|
2051 | else |
---|
2052 | { |
---|
2053 | sl2Act->next = new facet(2); |
---|
2054 | sl2Act = sl2Act->next; |
---|
2055 | sl2Act->setFacetNormal(f2Normal); |
---|
2056 | } |
---|
2057 | delete f2Normal; |
---|
2058 | codim2Act = codim2Act->next; |
---|
2059 | } |
---|
2060 | fAct = fAct->next; |
---|
2061 | delete fNormal; |
---|
2062 | }//if(fAct->isFlippable==TRUE) |
---|
2063 | else {fAct = fAct->next;} |
---|
2064 | }//End of copying facets into SLA |
---|
2065 | |
---|
2066 | SearchListAct = SearchListRoot; //Set to beginning of list |
---|
2067 | /*Make SearchList doubly linked*/ |
---|
2068 | while(SearchListAct!=NULL) |
---|
2069 | { |
---|
2070 | if(SearchListAct->next!=NULL) |
---|
2071 | { |
---|
2072 | SearchListAct->next->prev = SearchListAct; |
---|
2073 | } |
---|
2074 | SearchListAct = SearchListAct->next; |
---|
2075 | } |
---|
2076 | SearchListAct = SearchListRoot; //Set to beginning of List |
---|
2077 | |
---|
2078 | fAct = gcAct->facetPtr; |
---|
2079 | //NOTE Disabled until code works as expected. Reenabled 2.11.2009 |
---|
2080 | gcAct->writeConeToFile(*gcAct); |
---|
2081 | /*End of initialisation*/ |
---|
2082 | |
---|
2083 | fAct = SearchListAct; |
---|
2084 | /*2nd step |
---|
2085 | Choose a facet from fListPtr, flip it and forget the previous cone |
---|
2086 | We always choose the first facet from fListPtr as facet to be flipped |
---|
2087 | */ |
---|
2088 | while((SearchListAct!=NULL))// && counter<10) |
---|
2089 | {//NOTE See to it that the cone is only changed after ALL facets have been flipped! |
---|
2090 | fAct = SearchListAct; |
---|
2091 | |
---|
2092 | while(fAct!=NULL) |
---|
2093 | // while( (fAct->getUCN() == fAct->next->getUCN()) ) |
---|
2094 | { //Since SLA should only contain flippables there should be no need to check for that |
---|
2095 | gcAct->flip(gcAct->gcBasis,fAct); |
---|
2096 | //NOTE rCopy needed? |
---|
2097 | ring rTmp=rCopy(fAct->flipRing); |
---|
2098 | // ring rTmp=fAct->flipRing; //segfaults |
---|
2099 | rComplete(rTmp); |
---|
2100 | rChangeCurrRing(rTmp); |
---|
2101 | gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);//copy constructor! |
---|
2102 | /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing. |
---|
2103 | * Since we come at most once across a given facet from gcAct->facetPtr we can delete them. |
---|
2104 | * NOTE: Can this cause trouble with the destructor? |
---|
2105 | * Answer: Yes it bloody well does! |
---|
2106 | * However I'd like to delete this data here, since if we have a cone with many many facets it |
---|
2107 | * should pay to get rid of the flipGB as soon as possible. |
---|
2108 | * Destructor must be equipped with necessary checks. |
---|
2109 | */ |
---|
2110 | idDelete((ideal *)&fAct->flipGB); |
---|
2111 | rDelete(fAct->flipRing); |
---|
2112 | |
---|
2113 | gcTmp->getConeNormals(gcTmp->gcBasis, FALSE); |
---|
2114 | gcTmp->getCodim2Normals(*gcTmp); |
---|
2115 | gcTmp->normalize(); |
---|
2116 | //gcTmp->ddFacets should not be needed anymore, so |
---|
2117 | dd_FreeMatrix(gcTmp->ddFacets); |
---|
2118 | #ifdef gfan_DEBUG |
---|
2119 | // gcTmp->showFacets(1); |
---|
2120 | #endif |
---|
2121 | /*add facets to SLA here*/ |
---|
2122 | SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot); |
---|
2123 | if(gfanHeuristic==1) |
---|
2124 | { |
---|
2125 | gcTmp->writeConeToFile(*gcTmp); |
---|
2126 | //The for loop is no longer needed |
---|
2127 | // for(int ii=0;ii<IDELEMS(gcTmp->gcBasis);ii++) |
---|
2128 | // { |
---|
2129 | // pDelete(&gcTmp->gcBasis->m[ii]); |
---|
2130 | // } |
---|
2131 | idDelete((ideal*)&gcTmp->gcBasis);//Whonder why? |
---|
2132 | //If you use the following make sure it is uncommented in readConeFromFile |
---|
2133 | // rDelete(gcTmp->baseRing); |
---|
2134 | } |
---|
2135 | #ifdef gfan_DEBUG |
---|
2136 | // if(SearchListRoot!=NULL) |
---|
2137 | // gcTmp->showSLA(*SearchListRoot); |
---|
2138 | #endif |
---|
2139 | rChangeCurrRing(gcAct->baseRing); |
---|
2140 | rDelete(rTmp); |
---|
2141 | //doubly linked for easier removal |
---|
2142 | gcTmp->prev = gcPtr; |
---|
2143 | gcPtr->next=gcTmp; |
---|
2144 | gcPtr=gcPtr->next; |
---|
2145 | if(fAct->getUCN() == fAct->next->getUCN()) |
---|
2146 | { |
---|
2147 | fAct=fAct->next; |
---|
2148 | } |
---|
2149 | else |
---|
2150 | break; |
---|
2151 | // fAct=fAct->next; |
---|
2152 | }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) ); |
---|
2153 | //Search for cone with smallest UCN |
---|
2154 | gcNext = gcHead; |
---|
2155 | while(gcNext!=NULL && SearchListRoot!=NULL && gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && gcNext!=(gcone * const)0xfbfbfbfb) |
---|
2156 | { |
---|
2157 | if( gcNext->getUCN() == SearchListRoot->getUCN() ) |
---|
2158 | {//NOTE: Implement heuristic to be used! |
---|
2159 | if(gfanHeuristic==0) |
---|
2160 | { |
---|
2161 | gcAct = gcNext; |
---|
2162 | //Seems better to not to use rCopy here |
---|
2163 | // rAct=rCopy(gcAct->baseRing); |
---|
2164 | rAct=gcAct->baseRing; |
---|
2165 | rComplete(rAct); |
---|
2166 | rChangeCurrRing(rAct); |
---|
2167 | break; |
---|
2168 | } |
---|
2169 | else if(gfanHeuristic==1) |
---|
2170 | { |
---|
2171 | gcone *gcDel; |
---|
2172 | gcDel = gcAct; |
---|
2173 | gcAct = gcNext; |
---|
2174 | // rDelete(gcDel->baseRing); |
---|
2175 | //Read st00f from file |
---|
2176 | //implant the GB into gcAct |
---|
2177 | readConeFromFile(gcAct->getUCN(), gcAct); |
---|
2178 | // rAct=rCopy(gcAct->baseRing); |
---|
2179 | rAct=gcAct->baseRing; |
---|
2180 | rComplete(rAct); |
---|
2181 | rChangeCurrRing(rAct); |
---|
2182 | break; |
---|
2183 | } |
---|
2184 | } |
---|
2185 | gcNext = gcNext->next; |
---|
2186 | } |
---|
2187 | UCNcounter++; |
---|
2188 | //SearchListAct = SearchListAct->next; |
---|
2189 | //SearchListAct = fAct->next; |
---|
2190 | SearchListAct = SearchListRoot; |
---|
2191 | } |
---|
2192 | cout << endl << "Found " << counter << " cones - terminating" << endl; |
---|
2193 | }//void noRevS(gcone &gc) |
---|
2194 | |
---|
2195 | |
---|
2196 | /** \brief Make a set of rational vectors into integers |
---|
2197 | * |
---|
2198 | * We compute the lcm of the denominators and multiply with this to get integer values. |
---|
2199 | * \param dd_MatrixPtr,intvec |
---|
2200 | */ |
---|
2201 | inline void gcone::makeInt(const dd_MatrixPtr &M, const int line, intvec &n) |
---|
2202 | { |
---|
2203 | // mpz_t denom[this->numVars]; |
---|
2204 | mpz_t *denom = new mpz_t[this->numVars]; |
---|
2205 | for(int ii=0;ii<this->numVars;ii++) |
---|
2206 | { |
---|
2207 | mpz_init_set_str(denom[ii],"0",10); |
---|
2208 | } |
---|
2209 | |
---|
2210 | mpz_t kgV,tmp; |
---|
2211 | mpz_init(kgV); |
---|
2212 | mpz_init(tmp); |
---|
2213 | |
---|
2214 | for (int ii=0;ii<(M->colsize)-1;ii++) |
---|
2215 | { |
---|
2216 | mpz_t z; |
---|
2217 | mpz_init(z); |
---|
2218 | mpq_get_den(z,M->matrix[line-1][ii+1]); |
---|
2219 | mpz_set( denom[ii], z); |
---|
2220 | mpz_clear(z); |
---|
2221 | } |
---|
2222 | |
---|
2223 | /*Compute lcm of the denominators*/ |
---|
2224 | mpz_set(tmp,denom[0]); |
---|
2225 | for (int ii=0;ii<(M->colsize)-1;ii++) |
---|
2226 | { |
---|
2227 | mpz_lcm(kgV,tmp,denom[ii]); |
---|
2228 | mpz_set(tmp,kgV); |
---|
2229 | } |
---|
2230 | mpz_clear(tmp); |
---|
2231 | /*Multiply the nominators by kgV*/ |
---|
2232 | mpq_t qkgV,res; |
---|
2233 | mpq_init(qkgV); |
---|
2234 | mpq_set_str(qkgV,"1",10); |
---|
2235 | mpq_canonicalize(qkgV); |
---|
2236 | |
---|
2237 | mpq_init(res); |
---|
2238 | mpq_set_str(res,"1",10); |
---|
2239 | mpq_canonicalize(res); |
---|
2240 | |
---|
2241 | mpq_set_num(qkgV,kgV); |
---|
2242 | |
---|
2243 | // mpq_canonicalize(qkgV); |
---|
2244 | for (int ii=0;ii<(M->colsize)-1;ii++) |
---|
2245 | { |
---|
2246 | mpq_mul(res,qkgV,M->matrix[line-1][ii+1]); |
---|
2247 | n[ii]=(int)mpz_get_d(mpq_numref(res)); |
---|
2248 | } |
---|
2249 | delete [] denom; |
---|
2250 | mpz_clear(kgV); |
---|
2251 | mpq_clear(qkgV); mpq_clear(res); |
---|
2252 | |
---|
2253 | } |
---|
2254 | /** |
---|
2255 | * For all codim-2-facets we compute the gcd of the components of the facet normal and |
---|
2256 | * divide it out. Thus we get a normalized representation of each |
---|
2257 | * (codim-2)-facet normal, i.e. a primitive vector |
---|
2258 | */ |
---|
2259 | inline void gcone::normalize() |
---|
2260 | { |
---|
2261 | int *ggT = new int; |
---|
2262 | *ggT=1; |
---|
2263 | facet *fAct; |
---|
2264 | facet *codim2Act; |
---|
2265 | fAct = this->facetPtr; |
---|
2266 | codim2Act = fAct->codim2Ptr; |
---|
2267 | // intvec *n = new intvec(this->numVars); |
---|
2268 | |
---|
2269 | //while(codim2Act->next!=NULL) |
---|
2270 | while(fAct!=NULL) |
---|
2271 | { |
---|
2272 | while(codim2Act!=NULL) |
---|
2273 | { |
---|
2274 | intvec *n; |
---|
2275 | n=codim2Act->getFacetNormal(); |
---|
2276 | for(int ii=0;ii<this->numVars;ii++) |
---|
2277 | { |
---|
2278 | *ggT = intgcd((*ggT),(*n)[ii]); |
---|
2279 | } |
---|
2280 | for(int ii=0;ii<this->numVars;ii++) |
---|
2281 | { |
---|
2282 | (*n)[ii] = ((*n)[ii])/(*ggT); |
---|
2283 | } |
---|
2284 | codim2Act->setFacetNormal(n); |
---|
2285 | codim2Act = codim2Act->next; |
---|
2286 | delete n; |
---|
2287 | } |
---|
2288 | fAct = fAct->next; |
---|
2289 | } |
---|
2290 | delete ggT; |
---|
2291 | // delete n; |
---|
2292 | |
---|
2293 | } |
---|
2294 | /** \brief Enqueue new facets into the searchlist |
---|
2295 | * The searchlist (SLA for short) is implemented as a FIFO |
---|
2296 | * Checks are done as follows: |
---|
2297 | * 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it. |
---|
2298 | * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the |
---|
2299 | * facet from SLA and do not add fAct. |
---|
2300 | * It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we |
---|
2301 | * can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA. |
---|
2302 | * Takes ptr to search list root |
---|
2303 | * Returns a pointer to new first element of Searchlist |
---|
2304 | */ |
---|
2305 | facet * gcone::enqueueNewFacets(facet *f) |
---|
2306 | { |
---|
2307 | #ifdef gfanp |
---|
2308 | timeval start, end; |
---|
2309 | gettimeofday(&start, 0); |
---|
2310 | #endif |
---|
2311 | facet *slHead; |
---|
2312 | slHead = f; |
---|
2313 | facet *slAct; //called with f=SearchListRoot |
---|
2314 | slAct = f; |
---|
2315 | facet *slEnd; //Pointer to end of SLA |
---|
2316 | slEnd = f; |
---|
2317 | // facet *slEndStatic; //marks the end before a new facet is added |
---|
2318 | facet *fAct; |
---|
2319 | fAct = this->facetPtr; |
---|
2320 | facet *codim2Act; |
---|
2321 | codim2Act = this->facetPtr->codim2Ptr; |
---|
2322 | facet *sl2Act; |
---|
2323 | sl2Act = f->codim2Ptr; |
---|
2324 | /** Pointer to a facet that will be deleted */ |
---|
2325 | volatile facet *deleteMarker; |
---|
2326 | deleteMarker = NULL; |
---|
2327 | |
---|
2328 | /** Flag to indicate a facet that should be added to SLA*/ |
---|
2329 | // bool doNotAdd=FALSE; |
---|
2330 | /** \brief Flag to mark a facet that might be added |
---|
2331 | * The following case may occur: |
---|
2332 | * A facet fAct is found to be parallel but not equal to the current facet slAct from SLA. |
---|
2333 | * This does however not mean that there does not exist a facet behind the current slAct that is actually equal |
---|
2334 | * to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to |
---|
2335 | * FALSE and the according slAct is deleted. |
---|
2336 | * If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added |
---|
2337 | */ |
---|
2338 | // volatile bool maybe=FALSE; |
---|
2339 | /**A facet was removed, lengthOfSearchlist-- thus we must not rely on |
---|
2340 | * if(notParallelCtr==lengthOfSearchList) but rather |
---|
2341 | * if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) |
---|
2342 | */ |
---|
2343 | volatile bool removalOccured=FALSE; |
---|
2344 | int ctr=0; //encountered equalities in SLA |
---|
2345 | int notParallelCtr=0; |
---|
2346 | int lengthOfSearchList=1; |
---|
2347 | while(slEnd->next!=NULL) |
---|
2348 | { |
---|
2349 | slEnd=slEnd->next; |
---|
2350 | lengthOfSearchList++; |
---|
2351 | } |
---|
2352 | /*1st step: compare facetNormals*/ |
---|
2353 | // intvec *fNormal=NULL; |
---|
2354 | // intvec *slNormal=NULL; |
---|
2355 | |
---|
2356 | while(fAct!=NULL) |
---|
2357 | { |
---|
2358 | if(fAct->isFlippable==TRUE) |
---|
2359 | { |
---|
2360 | intvec *fNormal=NULL; |
---|
2361 | // intvec *slNormal=NULL; |
---|
2362 | fNormal=fAct->getFacetNormal(); |
---|
2363 | slAct = slHead; |
---|
2364 | notParallelCtr=0; |
---|
2365 | /*If slAct==NULL and fAct!=NULL |
---|
2366 | we just copy all remaining facets into SLA*/ |
---|
2367 | if(slAct==NULL) |
---|
2368 | { |
---|
2369 | facet *fCopy; |
---|
2370 | fCopy = fAct; |
---|
2371 | while(fCopy!=NULL) |
---|
2372 | { |
---|
2373 | if(slAct==NULL) |
---|
2374 | { |
---|
2375 | slAct = new facet(*fCopy); |
---|
2376 | slHead = slAct; |
---|
2377 | } |
---|
2378 | else |
---|
2379 | { |
---|
2380 | slAct->next = new facet(*fCopy); |
---|
2381 | slAct = slAct->next; |
---|
2382 | } |
---|
2383 | fCopy = fCopy->next; |
---|
2384 | } |
---|
2385 | break; |
---|
2386 | } |
---|
2387 | /*End of dumping into SLA*/ |
---|
2388 | while(slAct!=NULL) |
---|
2389 | //while(slAct!=slEndStatic->next) |
---|
2390 | { |
---|
2391 | intvec *slNormal=NULL; |
---|
2392 | removalOccured=FALSE; |
---|
2393 | slNormal = slAct->getFacetNormal(); |
---|
2394 | #ifdef gfan_DEBUG |
---|
2395 | cout << "Checking facet ("; |
---|
2396 | fNormal->show(1,1); |
---|
2397 | cout << ") against ("; |
---|
2398 | slNormal->show(1,1); |
---|
2399 | cout << ")" << endl; |
---|
2400 | #endif |
---|
2401 | if(areEqual(fAct,slAct)) |
---|
2402 | { |
---|
2403 | deleteMarker = slAct; |
---|
2404 | if(slAct==slHead) |
---|
2405 | { |
---|
2406 | slHead = slAct->next; |
---|
2407 | if(slHead!=NULL) |
---|
2408 | slHead->prev = NULL; |
---|
2409 | } |
---|
2410 | else if (slAct==slEnd) |
---|
2411 | { |
---|
2412 | slEnd=slEnd->prev; |
---|
2413 | slEnd->next = NULL; |
---|
2414 | } |
---|
2415 | else |
---|
2416 | { |
---|
2417 | slAct->prev->next = slAct->next; |
---|
2418 | slAct->next->prev = slAct->prev; |
---|
2419 | } |
---|
2420 | removalOccured=TRUE; |
---|
2421 | lengthOfSearchList--; |
---|
2422 | if(deleteMarker!=NULL) |
---|
2423 | { |
---|
2424 | // delete deleteMarker; |
---|
2425 | // deleteMarker=NULL; |
---|
2426 | } |
---|
2427 | #ifdef gfan_DEBUG |
---|
2428 | cout << "Removing ("; |
---|
2429 | fNormal->show(1,1); |
---|
2430 | cout << ") from list" << endl; |
---|
2431 | #endif |
---|
2432 | delete slNormal; |
---|
2433 | break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct |
---|
2434 | } |
---|
2435 | |
---|
2436 | slAct = slAct->next; |
---|
2437 | /* NOTE The following lines must not be here but rather called inside the if-block above. |
---|
2438 | Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now, |
---|
2439 | (not nice!) but since they are in seperate branches of the if-statement there should not |
---|
2440 | be a way it gets called twice thus ommiting one facet: |
---|
2441 | slAct = slAct->next;*/ |
---|
2442 | if(deleteMarker!=NULL) |
---|
2443 | { |
---|
2444 | // delete deleteMarker; |
---|
2445 | // deleteMarker=NULL; |
---|
2446 | } |
---|
2447 | delete slNormal; |
---|
2448 | //if slAct was marked as to be deleted, delete it here! |
---|
2449 | }//while(slAct!=NULL) |
---|
2450 | if(removalOccured==FALSE) |
---|
2451 | // if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) || (doNotAdd==FALSE) ) |
---|
2452 | //if( (notParallelCtr==lengthOfSearchList ) || doNotAdd==FALSE ) |
---|
2453 | { |
---|
2454 | #ifdef gfan_DEBUG |
---|
2455 | cout << "Adding facet ("; |
---|
2456 | fNormal->show(1,0); |
---|
2457 | cout << ") to SLA " << endl; |
---|
2458 | #endif |
---|
2459 | //Add fAct to SLA |
---|
2460 | facet *marker; |
---|
2461 | marker = slEnd; |
---|
2462 | facet *f2Act; |
---|
2463 | f2Act = fAct->codim2Ptr; |
---|
2464 | |
---|
2465 | slEnd->next = new facet(); |
---|
2466 | slEnd = slEnd->next; |
---|
2467 | facet *slEndCodim2Root; |
---|
2468 | facet *slEndCodim2Act; |
---|
2469 | slEndCodim2Root = slEnd->codim2Ptr; |
---|
2470 | slEndCodim2Act = slEndCodim2Root; |
---|
2471 | |
---|
2472 | slEnd->setUCN(this->getUCN()); |
---|
2473 | slEnd->isFlippable = TRUE; |
---|
2474 | slEnd->setFacetNormal(fNormal); |
---|
2475 | slEnd->prev = marker; |
---|
2476 | //Copy codim2-facets |
---|
2477 | // intvec *f2Normal=new intvec(this->numVars); |
---|
2478 | while(f2Act!=NULL) |
---|
2479 | { |
---|
2480 | intvec *f2Normal; |
---|
2481 | f2Normal=f2Act->getFacetNormal(); |
---|
2482 | if(slEndCodim2Root==NULL) |
---|
2483 | { |
---|
2484 | slEndCodim2Root = new facet(2); |
---|
2485 | slEnd->codim2Ptr = slEndCodim2Root; |
---|
2486 | slEndCodim2Root->setFacetNormal(f2Normal); |
---|
2487 | slEndCodim2Act = slEndCodim2Root; |
---|
2488 | } |
---|
2489 | else |
---|
2490 | { |
---|
2491 | slEndCodim2Act->next = new facet(2); |
---|
2492 | slEndCodim2Act = slEndCodim2Act->next; |
---|
2493 | slEndCodim2Act->setFacetNormal(f2Normal); |
---|
2494 | } |
---|
2495 | f2Act = f2Act->next; |
---|
2496 | delete f2Normal; |
---|
2497 | } |
---|
2498 | lengthOfSearchList++; |
---|
2499 | }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) || |
---|
2500 | fAct = fAct->next; |
---|
2501 | delete fNormal; |
---|
2502 | // delete slNormal; |
---|
2503 | }//if(fAct->isFlippable==TRUE) |
---|
2504 | else |
---|
2505 | { |
---|
2506 | fAct = fAct->next; |
---|
2507 | } |
---|
2508 | }//while(fAct!=NULL) |
---|
2509 | #ifdef gfanp |
---|
2510 | gettimeofday(&end, 0); |
---|
2511 | time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec)); |
---|
2512 | #endif |
---|
2513 | return slHead; |
---|
2514 | }//addC2N |
---|
2515 | |
---|
2516 | /** \brief Compute the gcd of two ints |
---|
2517 | */ |
---|
2518 | inline int gcone::intgcd(const int a, const int b) |
---|
2519 | { |
---|
2520 | int r, p=a, q=b; |
---|
2521 | if(p < 0) |
---|
2522 | p = -p; |
---|
2523 | if(q < 0) |
---|
2524 | q = -q; |
---|
2525 | while(q != 0) |
---|
2526 | { |
---|
2527 | r = p % q; |
---|
2528 | p = q; |
---|
2529 | q = r; |
---|
2530 | } |
---|
2531 | return p; |
---|
2532 | } |
---|
2533 | |
---|
2534 | /** \brief Construct a dd_MatrixPtr from a cone's list of facets |
---|
2535 | * NO LONGER USED |
---|
2536 | */ |
---|
2537 | inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc) |
---|
2538 | { |
---|
2539 | facet *fAct; |
---|
2540 | fAct = gc.facetPtr; |
---|
2541 | |
---|
2542 | dd_MatrixPtr M; |
---|
2543 | dd_rowrange ddrows; |
---|
2544 | dd_colrange ddcols; |
---|
2545 | ddcols=(this->numVars)+1; |
---|
2546 | ddrows=this->numFacets; |
---|
2547 | dd_NumberType numb = dd_Integer; |
---|
2548 | M=dd_CreateMatrix(ddrows,ddcols); |
---|
2549 | |
---|
2550 | int jj=0; |
---|
2551 | |
---|
2552 | while(fAct!=NULL) |
---|
2553 | { |
---|
2554 | intvec *comp; |
---|
2555 | comp = fAct->getFacetNormal(); |
---|
2556 | for(int ii=0;ii<this->numVars;ii++) |
---|
2557 | { |
---|
2558 | dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]); |
---|
2559 | } |
---|
2560 | jj++; |
---|
2561 | delete comp; |
---|
2562 | fAct=fAct->next; |
---|
2563 | } |
---|
2564 | return M; |
---|
2565 | } |
---|
2566 | |
---|
2567 | /** \brief Write information about a cone into a file on disk |
---|
2568 | * |
---|
2569 | * This methods writes the information needed for the "second" method into a file. |
---|
2570 | * The file's is divided in sections containing information on |
---|
2571 | * 1) the ring |
---|
2572 | * 2) the cone's Groebner Basis |
---|
2573 | * 3) the cone's facets |
---|
2574 | * Each line contains exactly one date |
---|
2575 | * Each section starts with its name in CAPITALS |
---|
2576 | */ |
---|
2577 | inline void gcone::writeConeToFile(const gcone &gc, bool usingIntPoints) |
---|
2578 | { |
---|
2579 | int UCN=gc.UCN; |
---|
2580 | stringstream ss; |
---|
2581 | ss << UCN; |
---|
2582 | string UCNstr = ss.str(); |
---|
2583 | |
---|
2584 | string prefix="/tmp/cone"; |
---|
2585 | // string prefix="./"; //crude hack -> run tests in separate directories |
---|
2586 | string suffix=".sg"; |
---|
2587 | string filename; |
---|
2588 | filename.append(prefix); |
---|
2589 | filename.append(UCNstr); |
---|
2590 | filename.append(suffix); |
---|
2591 | |
---|
2592 | |
---|
2593 | // int thisPID = getpid(); |
---|
2594 | // ss << thisPID; |
---|
2595 | // string strPID = ss.str(); |
---|
2596 | // string prefix="./"; |
---|
2597 | |
---|
2598 | ofstream gcOutputFile(filename.c_str()); |
---|
2599 | facet *fAct; |
---|
2600 | fAct = gc.facetPtr; |
---|
2601 | facet *f2Act; |
---|
2602 | f2Act=fAct->codim2Ptr; |
---|
2603 | |
---|
2604 | char *ringString = rString(currRing); |
---|
2605 | |
---|
2606 | if (!gcOutputFile) |
---|
2607 | { |
---|
2608 | cout << "Error opening file for writing in writeConeToFile" << endl; |
---|
2609 | } |
---|
2610 | else |
---|
2611 | { |
---|
2612 | gcOutputFile << "UCN" << endl; |
---|
2613 | gcOutputFile << this->UCN << endl; |
---|
2614 | gcOutputFile << "RING" << endl; |
---|
2615 | gcOutputFile << ringString << endl; |
---|
2616 | gcOutputFile << "GCBASISLENGTH" << endl; |
---|
2617 | gcOutputFile << IDELEMS(gc.gcBasis) << endl; |
---|
2618 | //Write this->gcBasis into file |
---|
2619 | gcOutputFile << "GCBASIS" << endl; |
---|
2620 | for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++) |
---|
2621 | { |
---|
2622 | gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl; |
---|
2623 | } |
---|
2624 | |
---|
2625 | gcOutputFile << "FACETS" << endl; |
---|
2626 | |
---|
2627 | while(fAct!=NULL) |
---|
2628 | { |
---|
2629 | intvec *iv; |
---|
2630 | iv=fAct->getFacetNormal(); |
---|
2631 | f2Act=fAct->codim2Ptr; |
---|
2632 | for (int ii=0;ii<iv->length();ii++) |
---|
2633 | { |
---|
2634 | if (ii<iv->length()-1) |
---|
2635 | { |
---|
2636 | gcOutputFile << (*iv)[ii] << ","; |
---|
2637 | } |
---|
2638 | else |
---|
2639 | { |
---|
2640 | gcOutputFile << (*iv)[ii] << " "; |
---|
2641 | } |
---|
2642 | } |
---|
2643 | delete iv; |
---|
2644 | while(f2Act!=NULL) |
---|
2645 | { |
---|
2646 | intvec *iv2; |
---|
2647 | iv2=f2Act->getFacetNormal(); |
---|
2648 | for(int jj=0;jj<iv2->length();jj++) |
---|
2649 | { |
---|
2650 | if (jj<iv2->length()-1) |
---|
2651 | { |
---|
2652 | gcOutputFile << (*iv2)[jj] << ","; |
---|
2653 | } |
---|
2654 | else |
---|
2655 | { |
---|
2656 | gcOutputFile << (*iv2)[jj] << " "; |
---|
2657 | } |
---|
2658 | } |
---|
2659 | delete iv2; |
---|
2660 | f2Act = f2Act->next; |
---|
2661 | } |
---|
2662 | gcOutputFile << endl; |
---|
2663 | fAct=fAct->next; |
---|
2664 | } |
---|
2665 | gcOutputFile.close(); |
---|
2666 | } |
---|
2667 | |
---|
2668 | }//writeConeToFile(gcone const &gc) |
---|
2669 | |
---|
2670 | /** \brief Reads a cone from a file identified by its number |
---|
2671 | */ |
---|
2672 | inline void gcone::readConeFromFile(int UCN, gcone *gc) |
---|
2673 | { |
---|
2674 | //int UCN=gc.UCN; |
---|
2675 | stringstream ss; |
---|
2676 | ss << UCN; |
---|
2677 | string UCNstr = ss.str(); |
---|
2678 | // string line; |
---|
2679 | // string strGcBasisLength; |
---|
2680 | // string strMonom, strCoeff, strCoeffNom, strCoeffDenom; |
---|
2681 | int gcBasisLength=0; |
---|
2682 | // int intCoeff=1; |
---|
2683 | // int intCoeffNom=1; //better (number) but that's not compatible with stringstream; |
---|
2684 | // int intCoeffDenom=1; |
---|
2685 | |
---|
2686 | // bool hasCoeffInQ = FALSE; //does the polynomial have rational coeff? |
---|
2687 | // bool hasNegCoeff = FALSE; //or a negative one? |
---|
2688 | size_t found; //used for find_first_(not)_of |
---|
2689 | |
---|
2690 | string prefix="/tmp/cone"; |
---|
2691 | string suffix=".sg"; |
---|
2692 | string filename; |
---|
2693 | filename.append(prefix); |
---|
2694 | filename.append(UCNstr); |
---|
2695 | filename.append(suffix); |
---|
2696 | |
---|
2697 | ifstream gcInputFile(filename.c_str(), ifstream::in); |
---|
2698 | |
---|
2699 | ring saveRing=currRing; |
---|
2700 | //Comment the following if you uncomment the if(line=="RING") part below |
---|
2701 | rChangeCurrRing(gc->baseRing); |
---|
2702 | |
---|
2703 | while( !gcInputFile.eof() ) |
---|
2704 | { |
---|
2705 | string line; |
---|
2706 | getline(gcInputFile,line); |
---|
2707 | // hasCoeffInQ = FALSE; |
---|
2708 | // hasNegCoeff = FALSE; |
---|
2709 | |
---|
2710 | if(line=="RING") |
---|
2711 | { |
---|
2712 | // getline(gcInputFile,line); |
---|
2713 | // found = line.find("a("); |
---|
2714 | // line.erase(0,found+2); |
---|
2715 | // string strweight; |
---|
2716 | // strweight=line.substr(0,line.find_first_of(")")); |
---|
2717 | // intvec *iv=new intvec(this->numVars); |
---|
2718 | // for(int ii=0;ii<this->numVars;ii++) |
---|
2719 | // { |
---|
2720 | // string weight; |
---|
2721 | // weight=line.substr(0,line.find_first_of(",)")); |
---|
2722 | // (*iv)[ii]=atoi(weight.c_str()); |
---|
2723 | // line.erase(0,line.find_first_of(",)")+1); |
---|
2724 | // } |
---|
2725 | // ring newRing; |
---|
2726 | // if(currRing->order[0]!=ringorder_a) |
---|
2727 | // { |
---|
2728 | // newRing=rCopyAndAddWeight(currRing,iv); |
---|
2729 | // } |
---|
2730 | // else |
---|
2731 | // { |
---|
2732 | // newRing=rCopy0(currRing); |
---|
2733 | // int length=this->numVars; |
---|
2734 | // int *A=(int *)omAlloc0(length*sizeof(int)); |
---|
2735 | // for(int jj=0;jj<length;jj++) |
---|
2736 | // { |
---|
2737 | // A[jj]=-(*iv)[jj]; |
---|
2738 | // } |
---|
2739 | // omFree(newRing->wvhdl[0]); |
---|
2740 | // newRing->wvhdl[0]=(int*)A; |
---|
2741 | // newRing->block1[0]=length; |
---|
2742 | // } |
---|
2743 | // delete iv; |
---|
2744 | // rComplete(newRing); |
---|
2745 | // gc->baseRing=rCopy(newRing); |
---|
2746 | // if(currRing!=gc->baseRing) |
---|
2747 | // rChangeCurrRing(gc->baseRing); |
---|
2748 | } |
---|
2749 | |
---|
2750 | if(line=="GCBASISLENGTH") |
---|
2751 | { |
---|
2752 | string strGcBasisLength; |
---|
2753 | getline(gcInputFile, line); |
---|
2754 | strGcBasisLength = line; |
---|
2755 | int size=atoi(strGcBasisLength.c_str()); |
---|
2756 | gcBasisLength=size; |
---|
2757 | gc->gcBasis=idInit(size,1); |
---|
2758 | } |
---|
2759 | if(line=="GCBASIS") |
---|
2760 | { |
---|
2761 | for(int jj=0;jj<gcBasisLength;jj++) |
---|
2762 | { |
---|
2763 | getline(gcInputFile,line); |
---|
2764 | //magically convert strings into polynomials |
---|
2765 | //polys.cc:p_Read |
---|
2766 | //check until first occurance of + or - |
---|
2767 | //data or c_str |
---|
2768 | poly strPoly=pInit();//Ought to be inside the while loop, but that will eat your memory |
---|
2769 | poly resPoly=pInit(); //The poly to be read in |
---|
2770 | while(!line.empty()) |
---|
2771 | { |
---|
2772 | // poly strPoly=pInit(); |
---|
2773 | number nCoeff=nInit(1); |
---|
2774 | number nCoeffNom=nInit(1); |
---|
2775 | number nCoeffDenom=nInit(1); |
---|
2776 | string strMonom, strCoeff, strCoeffNom, strCoeffDenom; |
---|
2777 | bool hasCoeffInQ = FALSE; //does the polynomial have rational coeff? |
---|
2778 | bool hasNegCoeff = FALSE; //or a negative one? |
---|
2779 | found = line.find_first_of("+-"); //get the first monomial |
---|
2780 | string tmp; |
---|
2781 | tmp=line[found]; |
---|
2782 | // if(found!=0 && (tmp.compare("-")==0) ) |
---|
2783 | // hasNegCoeff = TRUE; //We have a coeff < 0 |
---|
2784 | if(found==0) |
---|
2785 | { |
---|
2786 | if(tmp.compare("-")==0) |
---|
2787 | hasNegCoeff = TRUE; |
---|
2788 | line.erase(0,1); //remove leading + or - |
---|
2789 | found = line.find_first_of("+-"); //adjust found |
---|
2790 | } |
---|
2791 | strMonom = line.substr(0,found); |
---|
2792 | line.erase(0,found); |
---|
2793 | found = strMonom.find_first_of("/"); |
---|
2794 | if(found!=string::npos) //i.e. "/" exists in strMonom |
---|
2795 | { |
---|
2796 | hasCoeffInQ = TRUE; |
---|
2797 | strCoeffNom=strMonom.substr(0,found); |
---|
2798 | strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890",found+1)); |
---|
2799 | strMonom.erase(0,found); |
---|
2800 | strMonom.erase(0,strMonom.find_first_not_of("1234567890/")); |
---|
2801 | nRead(strCoeffNom.c_str(), &nCoeffNom); |
---|
2802 | nRead(strCoeffDenom.c_str(), &nCoeffDenom); |
---|
2803 | } |
---|
2804 | else |
---|
2805 | { |
---|
2806 | found = strMonom.find_first_not_of("1234567890"); |
---|
2807 | strCoeff = strMonom.substr(0,found); |
---|
2808 | if(!strCoeff.empty()) |
---|
2809 | { |
---|
2810 | nRead(strCoeff.c_str(),&nCoeff); |
---|
2811 | } |
---|
2812 | else |
---|
2813 | { |
---|
2814 | // intCoeff = 1; |
---|
2815 | nCoeff = nInit(1); |
---|
2816 | } |
---|
2817 | |
---|
2818 | } |
---|
2819 | const char* monom = strMonom.c_str(); |
---|
2820 | |
---|
2821 | p_Read(monom,strPoly,currRing); //strPoly:=monom |
---|
2822 | switch (hasCoeffInQ) |
---|
2823 | { |
---|
2824 | case TRUE: |
---|
2825 | if(hasNegCoeff) |
---|
2826 | nCoeffNom=nNeg(nCoeffNom); |
---|
2827 | pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom)); |
---|
2828 | break; |
---|
2829 | case FALSE: |
---|
2830 | if(hasNegCoeff) |
---|
2831 | nCoeff=nNeg(nCoeff); |
---|
2832 | if(!nIsOne(nCoeff)) |
---|
2833 | { |
---|
2834 | // if(hasNegCoeff) |
---|
2835 | // intCoeff *= -1; |
---|
2836 | // pSetCoeff(strPoly,(number) intCoeff); |
---|
2837 | pSetCoeff(strPoly, nCoeff ); |
---|
2838 | } |
---|
2839 | break; |
---|
2840 | |
---|
2841 | } |
---|
2842 | //pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1??? |
---|
2843 | if(pIsConstantComp(resPoly)) |
---|
2844 | resPoly=pCopy(strPoly); |
---|
2845 | else |
---|
2846 | resPoly=pAdd(resPoly,strPoly); |
---|
2847 | nDelete(&nCoeff); |
---|
2848 | nDelete(&nCoeffNom); |
---|
2849 | nDelete(&nCoeffDenom); |
---|
2850 | // pDelete(&strPoly); |
---|
2851 | }//while(!line.empty()) |
---|
2852 | gc->gcBasis->m[jj]=pCopy(resPoly); |
---|
2853 | pDelete(&resPoly); //reset |
---|
2854 | // pDelete(&strPoly); //NOTE Crashes |
---|
2855 | } |
---|
2856 | break; |
---|
2857 | }//if(line=="GCBASIS") |
---|
2858 | }//while(!gcInputFile.eof()) |
---|
2859 | gcInputFile.close(); |
---|
2860 | rChangeCurrRing(saveRing); |
---|
2861 | } |
---|
2862 | |
---|
2863 | ring gcone::getBaseRing() |
---|
2864 | { |
---|
2865 | return rCopy(this->baseRing); |
---|
2866 | } |
---|
2867 | /** \brief Gather the output |
---|
2868 | * List of lists |
---|
2869 | *\param *gc Pointer to gcone, preferably gcRoot ;-) |
---|
2870 | *\param n the number of cones |
---|
2871 | * |
---|
2872 | */ |
---|
2873 | lists lprepareResult(gcone *gc, int n) |
---|
2874 | { |
---|
2875 | gcone *gcAct; |
---|
2876 | gcAct = gc; |
---|
2877 | facet *fAct; |
---|
2878 | fAct = gc->facetPtr; |
---|
2879 | |
---|
2880 | lists res=(lists)omAllocBin(slists_bin); |
---|
2881 | res->Init(n); //initialize to store n cones |
---|
2882 | for(int ii=0;ii<n;ii++) |
---|
2883 | { |
---|
2884 | res->m[ii].rtyp=LIST_CMD; |
---|
2885 | lists l=(lists)omAllocBin(slists_bin); |
---|
2886 | l->Init(6); |
---|
2887 | l->m[0].rtyp=INT_CMD; |
---|
2888 | l->m[0].data=(void*)gcAct->getUCN(); |
---|
2889 | l->m[1].rtyp=IDEAL_CMD; |
---|
2890 | /*The following is necessary for leaves in the tree of cones |
---|
2891 | * Since we don't use them in the computation and gcBasis is |
---|
2892 | * set to (poly)NULL in noRevS we need to get this back here. |
---|
2893 | */ |
---|
2894 | // if(gcAct->gcBasis->m[0]==(poly)NULL) |
---|
2895 | if(gfanHeuristic==1) |
---|
2896 | gcAct->readConeFromFile(gcAct->getUCN(),gcAct); |
---|
2897 | // ring saveRing=currRing; |
---|
2898 | // ring tmpRing=gcAct->getBaseRing; |
---|
2899 | // rChangeCurrRing(tmpRing); |
---|
2900 | // l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing()); |
---|
2901 | l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing()); |
---|
2902 | // rChangeCurrRing(saveRing); |
---|
2903 | |
---|
2904 | l->m[2].rtyp=INTVEC_CMD; |
---|
2905 | intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr)); |
---|
2906 | l->m[2].data=(void*)ivCopy(&iv); |
---|
2907 | |
---|
2908 | l->m[3].rtyp=LIST_CMD; |
---|
2909 | lists lCodim2List = (lists)omAllocBin(slists_bin); |
---|
2910 | lCodim2List->Init(gcAct->numFacets); |
---|
2911 | fAct = gcAct->facetPtr;//fAct->codim2Ptr; |
---|
2912 | int jj=0; |
---|
2913 | while(fAct!=NULL && jj<gcAct->numFacets) |
---|
2914 | { |
---|
2915 | lCodim2List->m[jj].rtyp=INTVEC_CMD; |
---|
2916 | intvec ivC2=(gcAct->f2M(gcAct,fAct,2)); |
---|
2917 | lCodim2List->m[jj].data=(void*)ivCopy(&ivC2); |
---|
2918 | jj++; |
---|
2919 | fAct = fAct->next; |
---|
2920 | } |
---|
2921 | l->m[3].data=(void*)lCodim2List; |
---|
2922 | l->m[4].rtyp=RING_CMD; |
---|
2923 | l->m[4].data=(void*)(gcAct->getBaseRing()); |
---|
2924 | l->m[5].rtyp=INT_CMD; |
---|
2925 | l->m[5].data=(void*)gcAct->getPredUCN(); |
---|
2926 | res->m[ii].data=(void*)l; |
---|
2927 | gcAct = gcAct->next; |
---|
2928 | } |
---|
2929 | return res; |
---|
2930 | } |
---|
2931 | /** \brief Write facets of a cone into a matrix |
---|
2932 | * Takes a pointer to a facet as 2nd arg |
---|
2933 | * f should always point to gc->facetPtr |
---|
2934 | * param n is used to determine whether it operates in codim 1 or 2 |
---|
2935 | * |
---|
2936 | */ |
---|
2937 | inline intvec gcone::f2M(gcone *gc, facet *f, int n) |
---|
2938 | { |
---|
2939 | facet *fAct; |
---|
2940 | intvec *res;//=new intvec(this->numVars); |
---|
2941 | // int codim=n; |
---|
2942 | // int bound; |
---|
2943 | // if(f==gc->facetPtr) |
---|
2944 | if(n==1) |
---|
2945 | { |
---|
2946 | intvec *m1Res=new intvec(gc->numFacets,gc->numVars,0); |
---|
2947 | res = ivCopy(m1Res); |
---|
2948 | fAct = gc->facetPtr; |
---|
2949 | delete m1Res; |
---|
2950 | // bound = gc->numFacets*(this->numVars); |
---|
2951 | } |
---|
2952 | else |
---|
2953 | { |
---|
2954 | fAct = f->codim2Ptr; |
---|
2955 | intvec *m2Res = new intvec(f->numCodim2Facets,gc->numVars,0); |
---|
2956 | res = ivCopy(m2Res); |
---|
2957 | delete m2Res; |
---|
2958 | // bound = fAct->numCodim2Facets*(this->numVars); |
---|
2959 | |
---|
2960 | } |
---|
2961 | int ii=0; |
---|
2962 | while(fAct!=NULL )//&& ii < bound ) |
---|
2963 | { |
---|
2964 | intvec *fNormal; |
---|
2965 | fNormal = fAct->getFacetNormal(); |
---|
2966 | for(int jj=0;jj<this->numVars;jj++) |
---|
2967 | { |
---|
2968 | (*res)[ii]=(*fNormal)[jj]; |
---|
2969 | ii++; |
---|
2970 | } |
---|
2971 | fAct = fAct->next; |
---|
2972 | delete fNormal; |
---|
2973 | } |
---|
2974 | return *res; |
---|
2975 | } |
---|
2976 | |
---|
2977 | int gcone::counter=0; |
---|
2978 | int gfanHeuristic; |
---|
2979 | #ifdef gfanp |
---|
2980 | float gcone::time_getConeNormals; |
---|
2981 | float gcone::time_getCodim2Normals; |
---|
2982 | float gcone::time_flip; |
---|
2983 | float gcone::t_markings; |
---|
2984 | float gcone::t_dd; |
---|
2985 | float gcone::time_enqueue; |
---|
2986 | float gcone::time_computeInv; |
---|
2987 | float gcone::t_ddMC; |
---|
2988 | float gcone::t_mI; |
---|
2989 | float gcone::t_iP; |
---|
2990 | #endif |
---|
2991 | // ideal gfan(ideal inputIdeal, int h) |
---|
2992 | lists gfan(ideal inputIdeal, int h) |
---|
2993 | { |
---|
2994 | lists lResList; //this is the object we return |
---|
2995 | |
---|
2996 | if(rHasGlobalOrdering(currRing)) |
---|
2997 | { |
---|
2998 | int numvar = pVariables; |
---|
2999 | gfanHeuristic = h; |
---|
3000 | |
---|
3001 | enum searchMethod { |
---|
3002 | reverseSearch, |
---|
3003 | noRevS |
---|
3004 | }; |
---|
3005 | |
---|
3006 | searchMethod method; |
---|
3007 | method = noRevS; |
---|
3008 | |
---|
3009 | ring inputRing=currRing; // The ring the user entered |
---|
3010 | ring rootRing; // The ring associated to the target ordering |
---|
3011 | ideal res; |
---|
3012 | facet *fRoot; |
---|
3013 | dd_set_global_constants(); |
---|
3014 | if(method==noRevS) |
---|
3015 | { |
---|
3016 | gcone *gcRoot = new gcone(currRing,inputIdeal); |
---|
3017 | gcone *gcAct; |
---|
3018 | gcAct = gcRoot; |
---|
3019 | gcAct->numVars=pVariables; |
---|
3020 | gcAct->getGB(inputIdeal); |
---|
3021 | #ifndef NDEBUG |
---|
3022 | cout << "GB of input ideal is:" << endl; |
---|
3023 | // idShow(gcAct->gcBasis); |
---|
3024 | #endif |
---|
3025 | if(gcAct->isMonomial(gcAct->gcBasis)) |
---|
3026 | { |
---|
3027 | WerrorS("Monomial input - terminating"); |
---|
3028 | exit(-1); |
---|
3029 | gcAct->getConeNormals(gcAct->gcBasis); |
---|
3030 | lResList=lprepareResult(gcAct,1); |
---|
3031 | dd_free_global_constants; |
---|
3032 | //This is filthy |
---|
3033 | return lResList; |
---|
3034 | } |
---|
3035 | cout << endl; |
---|
3036 | gcAct->getConeNormals(gcAct->gcBasis); |
---|
3037 | gcAct->noRevS(*gcAct); |
---|
3038 | rChangeCurrRing(inputRing); |
---|
3039 | //res=gcAct->gcBasis; |
---|
3040 | //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS |
---|
3041 | res = inputIdeal; |
---|
3042 | lResList=lprepareResult(gcRoot,gcRoot->getCounter()); |
---|
3043 | /*Cleanup*/ |
---|
3044 | gcone *gcDel; |
---|
3045 | gcDel = gcRoot; |
---|
3046 | gcAct = gcRoot; |
---|
3047 | while(gcAct!=NULL) |
---|
3048 | { |
---|
3049 | gcDel = gcAct; |
---|
3050 | gcAct = gcAct->next; |
---|
3051 | delete gcDel; |
---|
3052 | } |
---|
3053 | } |
---|
3054 | dd_free_global_constants(); |
---|
3055 | }//rHasGlobalOrdering |
---|
3056 | else |
---|
3057 | { |
---|
3058 | WerrorS("Ring has non-global ordering - terminating"); |
---|
3059 | lResList->Init(1); |
---|
3060 | lResList->m[0].rtyp=INT_CMD; |
---|
3061 | int ires=0; |
---|
3062 | lResList->m[0].data=(void*)&ires; |
---|
3063 | } |
---|
3064 | //gcone::counter=0; |
---|
3065 | /*Return result*/ |
---|
3066 | #ifdef gfanp |
---|
3067 | cout << "t_getConeNormals:" << gcone::time_getConeNormals << endl; |
---|
3068 | cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl; |
---|
3069 | cout << " t_ddMC:" << gcone::t_ddMC << endl; |
---|
3070 | cout << " t_mI:" << gcone::t_mI << endl; |
---|
3071 | cout << " t_iP:" << gcone::t_iP << endl; |
---|
3072 | cout << "t_Flip:" << gcone::time_flip << endl; |
---|
3073 | cout << " t_markings:" << gcone::t_markings << endl; |
---|
3074 | cout << " t_dd:" << gcone::t_dd << endl; |
---|
3075 | cout << "t_computeInv:" << gcone::time_computeInv << endl; |
---|
3076 | cout << "t_enqueue:" << gcone::time_enqueue << endl; |
---|
3077 | #endif |
---|
3078 | return lResList; |
---|
3079 | } |
---|
3080 | |
---|
3081 | #endif |
---|