1 #include "wallideal.h"
2
3 #include "division.h"
4 #include "buchberger.h"
5 #include "timer.h"
6 #include "printer.h"
7 #include "lp.h"
8 #include "polyhedralcone.h"
9 #include "tropical2.h"
10 #include "linalg.h"
11 #include "log.h"
12
13 #include <set>
14 #include <iostream>
15
16 static Timer flipTimer("Flip",10);
17 static Timer flipTimer1("Flip1",10);
18 static Timer flipTimer2("Flip2",10);
19 static Timer flipTimer3("Flip3",10);
20 static Timer flipTimer4("Flip4",10);
21 static Timer flipTimer5("Flip5",10);
22 static Timer coneTimer("fajskfda",10);
23
wallPolynomial(Polynomial const & p,IntegerVector const & wallNormal)24 Polynomial wallPolynomial(Polynomial const &p, IntegerVector const &wallNormal)
25 {
26 Polynomial r(p.getRing());
27 IntegerVector markedExponent=p.getMarked().m.exponent;
28
29 for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
30 {
31 IntegerVector dif=markedExponent-i->first.exponent;
32 if(dependent(dif,wallNormal))
33 r+=Polynomial(Term(i->second,i->first));
34 }
35
36 r.mark(Monomial(r.getRing(),markedExponent));
37
38 return r;
39 }
40
wallPolynomial(Polynomial const & p,IntegerVectorList const & wallEqualities)41 static Polynomial wallPolynomial(Polynomial const &p, IntegerVectorList const &wallEqualities)
42 {
43 Polynomial r(p.getRing());
44 IntegerVector markedExponent=p.getMarked().m.exponent;
45
46 for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
47 {
48 IntegerVector dif=markedExponent-i->first.exponent;
49 bool dep=false;
50 for(IntegerVectorList::const_iterator j=wallEqualities.begin();j!=wallEqualities.end();j++)
51 {
52 if(dependent(dif,*j))
53 {
54 dep=true;
55 break;
56 }
57 }
58
59 if(dep || dif.isZero())
60 r+=Polynomial(Term(i->second,i->first));
61 }
62
63 r.mark(Monomial(p.getRing(),markedExponent));
64
65 return r;
66 }
67
wallIdeal(PolynomialSet const & groebnerBasis,IntegerVector const & wallNormal)68 PolynomialSet wallIdeal(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal)
69 {
70 PolynomialRing theRing=groebnerBasis.getRing();
71 PolynomialSet r(theRing);
72
73 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
74 {
75 r.push_back(wallPolynomial(*i,wallNormal));
76 }
77 return r;
78 }
79
lowerDimensionalWallIdeal(PolynomialSet const & groebnerBasis,IntegerVectorList const & wallEqualities)80 PolynomialSet lowerDimensionalWallIdeal(PolynomialSet const &groebnerBasis, IntegerVectorList const &wallEqualities)
81 {
82 PolynomialRing theRing=groebnerBasis.getRing();
83 PolynomialSet r(theRing);
84
85 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
86 {
87 r.push_back(wallPolynomial(*i,wallEqualities));
88 }
89 return r;
90 }
91
wallRemoveScaledInequalities(IntegerVectorList const & l)92 IntegerVectorList wallRemoveScaledInequalities(IntegerVectorList const &l)
93 {
94 IntegerVectorList ret;
95
96 /*for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
97 {
98 bool found=false;
99
100 assert(!i->isZero());
101
102 for(IntegerVectorList::const_iterator k=ret.begin();k!=ret.end();k++)
103 if(dependent(*i,*k)&&dotLong(*i,*k)>0)found=true;
104
105 if(!found)ret.push_back(*i);
106 }
107 */
108 set<IntegerVector> temp;
109 // log0 fprintf(Stderr,"C\n");
110 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)temp.insert(normalized(*i));
111 // log0 fprintf(Stderr,"D\n");
112 for(set<IntegerVector>::const_iterator i=temp.begin();i!=temp.end();i++)ret.push_back(*i);
113
114 return ret;
115 }
116
117
algebraicTest(IntegerVectorList const & l,PolynomialSet const & groebnerBasis)118 IntegerVectorList algebraicTest(IntegerVectorList const &l, PolynomialSet const &groebnerBasis) // TO DO: FIGURE OUT IF THIS TEST WORKS IN THE NON-HOMOGENEOUS CASE
119 {
120 PolynomialRing theRing=groebnerBasis.getRing();
121 LexicographicTermOrder T;
122 PolynomialSet LT(theRing);
123
124 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
125 {
126 LT.push_back(Polynomial(i->getMarked()));
127 }
128
129
130 //fprintf(Stderr,"In Size:%i\n",l.size());
131 IntegerVectorList ret;
132 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
133 {
134 bool accept=true;
135 PolynomialSet g2=wallIdeal(groebnerBasis,*i);
136 for(PolynomialSet::const_iterator j=g2.begin();j!=g2.end() && accept;j++)
137 for(PolynomialSet::const_iterator k=j;k!=g2.end();k++)
138 {
139 // fprintf(Stderr,"test\n");
140 if((!j->isMonomial()) || (!k->isMonomial()))
141 if(!relativelyPrime(j->getMarked().m.exponent,k->getMarked().m.exponent))
142 {
143 Polynomial s=sPolynomial(*j, *k);
144 if(!s.isZero())
145 {
146 bool witness=true;
147 for(TermMap::const_iterator i=s.terms.begin();i!=s.terms.end();i++)//TODO: prove that it suffices to check just the leading term of the S-polynomial and adjust the code accordingly.
148 {
149 bool inideal=false;
150 for(PolynomialSet::const_iterator j=LT.begin();j!=LT.end();j++)
151 if(j->getMarked().m.exponent.divides(i->first.exponent))
152 {
153 inideal=true;
154 break;
155 }
156 if(inideal)
157 {
158 witness=false;
159 break;
160 }
161 }
162 if(witness)
163 {
164 accept=false;
165 break;
166 }
167
168 /*
169 s.mark(T); // with respect to some termorder
170 s.scaleMarkedCoefficientToOne();
171
172 if(1)
173 {
174 Polynomial t=division(s,LT,T);
175 if((t-s).isZero())
176 {
177 accept=false;
178 break;
179 }
180 }
181 else
182 {
183 s=division(s,g2,T);
184 if(!s.isZero())
185 {
186 accept=false;
187 break;
188 }
189 }*/
190 }
191 }
192
193 }
194 if(accept)ret.push_back(*i);
195 }
196 // fprintf(Stderr,"Out Size:%i\n",ret.size());
197 return ret;
198 }
199
200
exponentDifferences(PolynomialSet const & groebnerBasis)201 IntegerVectorList exponentDifferences(PolynomialSet const &groebnerBasis)
202 {
203 IntegerVectorList ret;
204 set<IntegerVector> temp;
205 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
206 {
207 IntegerVector markedExponent=i->getMarked().m.exponent;
208
209 for(TermMap::const_iterator j=i->terms.begin();j!=i->terms.end();j++)
210 {
211 IntegerVector dif=markedExponent-j->first.exponent;
212
213 if(!dif.isZero())
214 {
215 /* bool found=false; //These four lines are now done in wallRemoveScaledInequalities()
216
217 for(IntegerVectorList::const_iterator k=ret.begin();k!=ret.end();k++)
218 if(dependent(dif,*k))found=true;
219
220 if(!found)
221 */
222 temp.insert(dif);
223 }
224 }
225 }
226 for(set<IntegerVector>::const_iterator i=temp.begin();i!=temp.end();i++)ret.push_back(*i);
227 return ret;
228 }
229
wallInequalities(PolynomialSet const & groebnerBasis)230 IntegerVectorList wallInequalities(PolynomialSet const &groebnerBasis)
231 {
232 return wallRemoveScaledInequalities(exponentDifferences(groebnerBasis));
233 // return algebraicTest(wallRemoveScaledInequalities(ret),groebnerBasis);
234 }
235
236
flip(PolynomialSet const & groebnerBasis,IntegerVector const & wallNormal,TermOrder * autoReduceHint)237 PolynomialSet flip(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal, TermOrder *autoReduceHint)
238 {
239 PolynomialRing theRing=groebnerBasis.getRing();
240 // fprintf(Stderr,"ENTERING flip\n");
241 // fprintf(Stderr,"flip - start\n");
242 // AsciiPrinter(Stderr).printPolynomialSet(groebnerBasis);
243 // AsciiPrinter(Stderr).printVector(wallNormal);
244
245 TimerScope ts(&flipTimer);
246 //Subroutine 3.7 in [Sturmfels]
247
248 // Step 1
249 // fprintf(Stderr,"flip - step1\n");
250 flipTimer1.on();
251 PolynomialSet wall=wallIdeal(groebnerBasis,wallNormal);
252 wall.markAndScale(WeightTermOrder(wallNormal));// This marking will be used later on when we lift
253 // fprintf(Stderr,"Changed order:\n");
254 // AsciiPrinter(Stderr).printPolynomialSet(wall);
255 flipTimer1.off();
256
257 // Step 2
258 // fprintf(Stderr,"flip - step2\n");
259 // AsciiPrinter(Stderr).printPolynomialSet(wall);
260 flipTimer2.on();
261 PolynomialSet oldWall=wall;
262 WeightTermOrder flipOrder(-wallNormal);
263 buchberger(&wall,flipOrder);
264 flipTimer2.off();
265
266 // Step 3
267 // fprintf(Stderr,"flip - step3\n");
268 flipTimer3.on();
269 PolynomialSet newBasis(theRing);
270 flipTimer3.off();
271
272 // Step 4 Lift
273 // fprintf(Stderr,"flip - lifting\n");
274 // fprintf(Stderr,"flip - step4\n");
275
276 {
277 // liftBasis() could be used for this!!!!
278
279 TimerScope ts(&flipTimer4);
280 for(PolynomialSet::const_iterator j=wall.begin();j!=wall.end();j++)
281 {
282 newBasis.push_back(divisionLift(*j, oldWall, groebnerBasis, LexicographicTermOrder()));
283 /*{
284 // The following should also work:
285 Polynomial q=*j-division(*j,groebnerBasis,LexicographicTermOrder());
286 assert(!q.isZero());
287
288 newBasis.push_back(q);
289 }*/
290 }
291 }
292
293 // Step 5 Autoreduce
294 // fprintf(Stderr,"flip - autoreduce\n");
295 // AsciiPrinter(Stderr).printPolynomialSet(wall);
296 // AsciiPrinter(Stderr).printPolynomialSet(newBasis);
297 // fprintf(Stderr,"flip - step5\n");
298 {
299 TimerScope ts(&flipTimer5);
300 PolynomialSet::const_iterator k=wall.begin();
301 for(PolynomialSet::iterator j=newBasis.begin();j!=newBasis.end();j++)
302 {
303 j->mark(k->getMarked().m);
304 k++;
305 }
306
307 // fprintf(Stderr,"Marked order:\n");
308 // AsciiPrinter(Stderr).printPolynomialSet(wall);
309 // AsciiPrinter(Stderr).printPolynomialSet(newBasis);
310 // fprintf(Stderr,"Not reduced lifted basis:\n");
311 // AsciiPrinter(Stderr).printPolynomialSet(newBasis);
312
313 static int t;
314 t++;
315 t&=0;
316 if(t==0)
317 {
318 // fprintf(Stderr,"autoreducing ..");
319 if(autoReduceHint)
320 autoReduce(&newBasis,*autoReduceHint);
321 else
322 autoReduce(&newBasis,LexicographicTermOrder());
323 // fprintf(Stderr,".. done\n");
324 }
325
326 //autoReduce(&newBasis,StandardGradedLexicographicTermOrder());
327 }
328
329 // fprintf(Stderr,"flip - done\n");
330 // fprintf(Stderr,"LEAVING flip\n");
331
332 //fprintf(Stderr,"%i",newBasis.size());
333
334 return newBasis;
335 }
336
337
wallContainsPositiveVector(IntegerVector const & wallNormal)338 bool wallContainsPositiveVector(IntegerVector const &wallNormal)
339 {
340 //This is not right I think
341 int n=wallNormal.size();
342 for(int i=0;i<n;i++)if(wallNormal[i]<0)return true;
343
344 return false;
345 }
346
wallAddCoordinateWalls(IntegerVectorList & normals)347 void wallAddCoordinateWalls(IntegerVectorList &normals)
348 {
349 assert(!normals.empty());
350 int n=normals.begin()->size();
351 for(int i=0;i<n;i++)normals.push_back(IntegerVector::standardVector(n,i));
352 }
353
354
isIdealHomogeneous(PolynomialSet const & groebnerBasis)355 bool isIdealHomogeneous(PolynomialSet const &groebnerBasis)
356 {
357 int n=groebnerBasis.numberOfVariablesInRing();
358 IntegerVectorList a;
359 PolyhedralCone p=intersection(PolyhedralCone(a,wallInequalities(groebnerBasis),n),PolyhedralCone::positiveOrthant(n));
360
361 return p.containsPositiveVector();
362 }
363
364
365 /* This routine is a preprocessing step for redudancy removal.
366 The routine normalizes the list of vectors in gcd sense.
367 It removes duplicates.
368 It removes direction that are sums of other directions.
369 The input must satisfy:
370 Input must be pointed, meaning that there must exist a codimension one subspace with all the input vectors strictly on one side.
371 It particular, there is no zero-vector in the list (It would be easy to change the routine to handle this case)
372 These requirements guarantee that *i and *j are not removed in the line with the comment.
373 */
374 /*
375 There are two versions of this routine.
376 */
377 //400 -> 80
normalizedWithSumsAndDuplicatesRemoved2(IntegerVectorList const & a)378 IntegerVectorList normalizedWithSumsAndDuplicatesRemoved2(IntegerVectorList const &a)
379 {
380 IntegerVectorList ret;
381 set<IntegerVector> b;
382
383 for(IntegerVectorList::const_iterator i=a.begin();i!=a.end();i++)
384 {
385 assert(!(i->isZero()));
386 b.insert(normalized(*i));
387 }
388
389 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
390 for(set<IntegerVector>::const_iterator j=i;j!=b.end();j++)
391 if(i!=j)b.erase(normalized(*i+*j));//this can never remove *i or *j
392
393 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
394 ret.push_back(*i);
395
396 return ret;
397 }
398
399
400 class MyHashMap
401 {
402 typedef vector<set<IntegerVector> > Container;
403 Container container;
404 int tableSize;
405 public:
406 class iterator
407 {
408 class MyHashMap &hashMap;
409 int index; // having index==-1 means that we are before/after the elements.
410 set<IntegerVector>::iterator i;
411 public:
operator ++()412 bool operator++()
413 {
414 if(index==-1)goto search;
415 i++;
416 while(i==hashMap.container[index].end())
417 {
418 search:
419 index++;
420 if(index>=hashMap.tableSize){
421 index=-1;
422 return false;
423 }
424 i=hashMap.container[index].begin();
425 }
426 return true;
427 }
operator *() const428 IntegerVector const & operator*()const
429 {
430 return *i;
431 }
operator *()432 IntegerVector operator*()
433 {
434 return *i;
435 }
436 /* bool operator()()
437 {
438 return index!=-1;
439 }*/
iterator(MyHashMap & hashMap_)440 iterator(MyHashMap &hashMap_):
441 hashMap(hashMap_)
442 {
443 index=-1;
444 }
445 };
function(const IntegerVector & v)446 unsigned int function(const IntegerVector &v)
447 {
448 unsigned int ret=0;
449 int n=v.size();
450 for(int i=0;i<n;i++)
451 ret=(ret<<3)+(ret>>29)+v.UNCHECKEDACCESS(i);
452 return ret%tableSize;
453 }
MyHashMap(int tableSize_)454 MyHashMap(int tableSize_):
455 container(tableSize_),
456 tableSize(tableSize_)
457 {
458 assert(tableSize_>0);
459 }
insert(const IntegerVector & v)460 void insert(const IntegerVector &v)
461 {
462 container[function(v)].insert(v);
463 }
erase(IntegerVector const & v)464 void erase(IntegerVector const &v)
465 {
466 container[function(v)].erase(v);
467 }
begin()468 iterator begin()
469 {
470 iterator ret(*this);
471 ++ret;
472 return ret;
473 }
size()474 int size()
475 {
476 iterator i=begin();
477 int ret=0;
478 do{ret++;}while(++i);
479 return ret;
480 }
print()481 void print()
482 {
483 for(int i=0;i<container.size();i++)
484 {
485 debug.printInteger(i);
486 debug<<":\n";
487 for(set<IntegerVector>::const_iterator j=container[i].begin();j!=container[i].end();j++)debug<<*j;
488 }
489 }
490 };
491
492
normalizedWithSumsAndDuplicatesRemoved(IntegerVectorList const & a)493 IntegerVectorList normalizedWithSumsAndDuplicatesRemoved(IntegerVectorList const &a)
494 {
495 if(a.empty())return a;
496 int n=a.front().size();
497 IntegerVector temp1(n);
498 IntegerVector temp2(n);
499 IntegerVectorList ret;
500 MyHashMap b(a.size());
501
502 // log0 debug<<"a.size"<<(int)a.size()<<"\n";
503 // log0 debug<<"b.size"<<(int)b.size()<<"\n";
504
505 // log0 cerr << "N:"<<a.size();
506 for(IntegerVectorList::const_iterator i=a.begin();i!=a.end();i++)
507 {
508 assert(!(i->isZero()));
509 b.insert(normalized(*i));
510 }
511 // b.print();
512 //debug<<"TEST";
513 // log0 cerr << "N:"<<b.size();
514
515 {
516 MyHashMap::iterator i=b.begin();
517
518 do
519 {
520 MyHashMap::iterator j=i;
521 while(++j)
522 {
523 // b.erase(normalized(*i+*j));//this can never remove *i or *j
524 IntegerVector const &I=*i;
525 IntegerVector const &J=*j;
526 for(int k=0;k<n;k++)temp1[k]=I.UNCHECKEDACCESS(k)+J.UNCHECKEDACCESS(k);
527 normalizedLowLevel(temp1,temp2);
528 b.erase(temp2);//this can never remove *i or *j
529 }
530 }
531 while(++i);
532 // log0 cerr << "N:"<<b.size();
533 }
534 vector<IntegerVector> original;
535 {
536 MyHashMap::iterator i=b.begin();
537
538 do
539 {
540 original.push_back(*i);
541 }
542 while(++i);
543 }
544
545 for(vector<IntegerVector>::const_iterator i=original.begin();i!=original.end();i++)
546 for(list<IntegerVector>::const_iterator j=a.begin();j!=a.end();j++)
547 /* for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
548 for(set<IntegerVector>::const_iterator j=i;j!=b.end();j++)*/
549 // if(*i!=*j)b.erase(normalized(*i+*j));//this can never remove *i or *j
550 if(!dependent(*i,*j))
551 {
552 // b.erase(normalized(*i+*j));//this can never remove *i or *j
553 IntegerVector const &I=*i;
554 IntegerVector const &J=*j;
555 for(int k=0;k<n;k++)temp1[k]=I.UNCHECKEDACCESS(k)+J.UNCHECKEDACCESS(k);
556 normalizedLowLevel(temp1,temp2);
557 b.erase(temp2);//this can never remove *i or *j
558 }
559
560 // log0 cerr << "N:"<<b.size();
561
562 // log0 debug<<"b.size"<<(int)b.size()<<"\n";
563 {
564 MyHashMap::iterator i=b.begin();
565
566 do
567 {
568 IntegerVector temp=*i;
569 ret.push_back(temp);
570 }
571 while(++i);
572 }
573 return ret;
574 }
575
576
577
578
579 // 400 -> 20
normalizedWithSumsAndDuplicatesRemovedNOHASH(IntegerVectorList const & a)580 IntegerVectorList normalizedWithSumsAndDuplicatesRemovedNOHASH(IntegerVectorList const &a)
581 {
582 /* IntegerVectorList a;
583 {
584 set<IntegerVector> b;
585 for(IntegerVectorList::const_iterator i=A.begin();i!=A.end();i++)b.insert(*i);
586 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)a.push_back(*i);
587 }*/
588
589 IntegerVectorList ret;
590 set<IntegerVector> b;
591
592 log0 cerr << "N:"<<a.size();
593 for(IntegerVectorList::const_iterator i=a.begin();i!=a.end();i++)
594 {
595 assert(!(i->isZero()));
596 b.insert(normalized(*i));
597 }
598
599 log0 cerr << "N:"<<b.size();
600
601 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
602 for(set<IntegerVector>::const_iterator j=i;j!=b.end();j++)
603 if(i!=j)b.erase(normalized(*i+*j));//this can never remove *i or *j
604
605 log0 cerr << "N:"<<b.size();
606
607 vector<IntegerVector> original;
608 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
609 original.push_back(*i);
610
611 for(vector<IntegerVector>::const_iterator i=original.begin();i!=original.end();i++)
612 for(list<IntegerVector>::const_iterator j=a.begin();j!=a.end();j++)
613 /* for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
614 for(set<IntegerVector>::const_iterator j=i;j!=b.end();j++)*/
615 // if(*i!=*j)b.erase(normalized(*i+*j));//this can never remove *i or *j
616 if(!dependent(*i,*j))b.erase(normalized(*i+*j));//this can never remove *i or *j
617
618 log0 cerr << "N:"<<b.size();
619
620 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
621 ret.push_back(*i);
622
623 return ret;
624 }
625
626
627 //400 -> 16
normalizedWithSumsAndDuplicatesRemoved3(IntegerVectorList const & a)628 IntegerVectorList normalizedWithSumsAndDuplicatesRemoved3(IntegerVectorList const &a)
629 {
630 IntegerVectorList ret;
631 set<IntegerVector> b;
632
633 for(IntegerVectorList::const_iterator i=a.begin();i!=a.end();i++)
634 {
635 assert(!(i->isZero()));
636 b.insert(normalized(*i));
637 }
638
639 vector<IntegerVector> original;
640 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
641 original.push_back(*i);
642
643 for(vector<IntegerVector>::const_iterator i=original.begin();i!=original.end();i++)
644 for(vector<IntegerVector>::const_iterator j=i;j!=original.end();j++)
645 /* for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
646 for(set<IntegerVector>::const_iterator j=i;j!=b.end();j++)*/
647 if(i!=j)b.erase(normalized(*i+*j));//this can never remove *i or *j
648
649 for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
650 ret.push_back(*i);
651
652 return ret;
653 }
654
655
wallFlipableNormals(PolynomialSet const & groebnerBasis,bool isKnownToBeHomogeneous)656 IntegerVectorList wallFlipableNormals(PolynomialSet const &groebnerBasis, bool isKnownToBeHomogeneous)
657 {
658 // New optimised version using PolyhedralCone
659 int n=groebnerBasis.numberOfVariablesInRing();
660 IntegerVectorList a;
661 //AsciiPrinter(Stderr).printVectorList(wallInequalities(groebnerBasis));
662 //PolyhedralCone p=intersection(PolyhedralCone(wallInequalities(groebnerBasis),a),PolyhedralCone::positiveOrthant(n));
663
664 IntegerVectorList normals=algebraicTest(wallInequalities(groebnerBasis),groebnerBasis);
665 coneTimer.on();
666 // PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
667
668 /* AsciiPrinter(Stderr).printVectorList(normals);
669 AsciiPrinter(Stderr).printInteger(normals.size());
670 */
671 normals=normalizedWithSumsAndDuplicatesRemoved(normals);
672
673 // AsciiPrinter(Stderr).printVectorList(normals);
674 // AsciiPrinter(Stderr).printInteger(normals.size());
675
676 PolyhedralCone p=PolyhedralCone(normals,a,n);
677 // fprintf(Stderr,"Finding facets\n");
678 p.findFacets();
679 // fprintf(Stderr,"Done Finding facets\n");
680 coneTimer.off();
681
682 IntegerVectorList ilist=p.getHalfSpaces();
683 IntegerVectorList ret;
684 for(IntegerVectorList::iterator i=ilist.begin();i!=ilist.end();i++)
685 {
686 bool facetOK=true;
687 if(!isKnownToBeHomogeneous)
688 {
689 *i=-1 * (*i);
690 PolyhedralCone temp=intersection(PolyhedralCone(ilist,a),PolyhedralCone::positiveOrthant(n));
691 facetOK=(temp.dimension()==n);
692 *i=-1 * (*i);
693 }
694
695 if(facetOK)
696 ret.push_back(*i);
697 }
698
699 // New version using PolyhedralCone
700 /* int n=groebnerBasis.numberOfVariablesInRing();
701 IntegerVectorList a;
702 // AsciiPrinter(Stderr).printVectorList(wallInequalities(groebnerBasis));
703 PolyhedralCone p=intersection(PolyhedralCone(wallInequalities(groebnerBasis),a),PolyhedralCone::positiveOrthant(n));
704 // PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
705 p.findFacets();
706 IntegerVectorList ilist=p.getHalfSpaces();
707 IntegerVectorList ret;
708 for(IntegerVectorList::const_iterator i=ilist.begin();i!=ilist.end();i++)
709 if(wallContainsPositiveVector(*i))ret.push_back(*i);
710 */
711
712 /* Old code not using PolyhedralCone
713 IntegerVectorList ilist=wallInequalities(groebnerBasis);
714 wallAddCoordinateWalls(ilist);
715 IntegerVectorList ret;
716 for(IntegerVectorList::const_iterator i=ilist.begin();i!=ilist.end();i++)
717 if(isFacet(ilist,i)&&wallContainsPositiveVector(*i))
718 ret.push_back(*i);
719 */
720 return ret;
721 }
722
723
flipMinkowski(Polynomial p,IntegerVector const & wallNormal)724 Polynomial flipMinkowski(Polynomial p, IntegerVector const &wallNormal)
725 {
726 IntegerVector best=p.getMarked().m.exponent;
727
728 for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
729 {
730 IntegerVector e=i->first.exponent;
731 IntegerVector diff=e-best;
732 if(dependent(diff,wallNormal)&&dot(diff,wallNormal)<0)best=e;
733 }
734 p.mark(Monomial(p.getRing(),best));
735
736 return p;
737 }
738
739
flipMinkowski(PolynomialSet const & groebnerBasis,IntegerVector const & wallNormal)740 PolynomialSet flipMinkowski(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal)
741 {
742 PolynomialSet r(groebnerBasis.getRing());
743
744 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
745 r.push_back(flipMinkowski(*i,wallNormal));
746
747 return r;
748 }
749
750
homogeneitySpace(PolynomialSet const & reducedGroebnerBasis)751 PolyhedralCone homogeneitySpace(PolynomialSet const &reducedGroebnerBasis)
752 {
753 IntegerVectorList l=wallInequalities(reducedGroebnerBasis);
754 IntegerVectorList a;
755 PolyhedralCone c(a,l,reducedGroebnerBasis.getRing().getNumberOfVariables());
756 c.findImpliedEquations();
757 // c.canonicalize();
758 return c;
759 }
760
groebnerCone(PolynomialSet const & reducedGroebnerBasis,bool useAlgebraicTest)761 PolyhedralCone groebnerCone(PolynomialSet const &reducedGroebnerBasis, bool useAlgebraicTest)
762 {
763 int n=reducedGroebnerBasis.getRing().getNumberOfVariables();
764 IntegerVectorList l=wallInequalities(reducedGroebnerBasis);
765 if(useAlgebraicTest)l=algebraicTest(l,reducedGroebnerBasis);
766 l=normalizedWithSumsAndDuplicatesRemoved(l);
767 IntegerVectorList a;
768 PolyhedralCone c(l,a,n);
769 c.canonicalize();
770 return c;
771 }
772
dimensionOfHomogeneitySpace(PolynomialSet const & reducedGroebnerBasis)773 int dimensionOfHomogeneitySpace(PolynomialSet const &reducedGroebnerBasis)
774 {
775 return homogeneitySpace(reducedGroebnerBasis).dimension();
776 }
777
778
liftBasis(PolynomialSet const & toBeLifted,PolynomialSet const & originalBasisForFullIdeal)779 PolynomialSet liftBasis(PolynomialSet const &toBeLifted, PolynomialSet const &originalBasisForFullIdeal)
780 {
781 PolynomialRing theRing=toBeLifted.getRing();
782 assert(toBeLifted.isValid());
783 assert(originalBasisForFullIdeal.isValid());
784
785 PolynomialSet newBasis(theRing);
786
787 // fprintf(Stderr,"LIFTING:");
788 // AsciiPrinter(Stderr).printPolynomialSet(toBeLifted);
789
790 for(PolynomialSet::const_iterator j=toBeLifted.begin();j!=toBeLifted.end();j++)
791 {
792 assert(!j->isZero());
793 //AsciiPrinter(Stderr).printVector(j->getMarked().m.exponent);
794
795 Polynomial q=*j-division(*j,originalBasisForFullIdeal,LexicographicTermOrder());
796 assert(!q.isZero());
797 q.mark(j->getMarked().m);
798 newBasis.push_back(q);
799 }
800 autoReduce(&newBasis,LexicographicTermOrder());
801 // fprintf(Stderr,"TO:");
802 // AsciiPrinter(Stderr).printPolynomialSet(newBasis);
803
804 return newBasis;
805 }
806
807
isMarkingConsistent(PolynomialSet const & g)808 bool isMarkingConsistent(PolynomialSet const &g)
809 {
810 IntegerVectorList empty;
811 PolyhedralCone c(wallInequalities(g),empty,g.getRing().getNumberOfVariables());
812 c=intersection(c,PolyhedralCone::positiveOrthant(c.ambientDimension()));
813 log1 AsciiPrinter(Stderr).printPolyhedralCone(c);
814 return c.dimension()==c.ambientDimension();
815 }
816
817
818 /*
819 Heuristic for checking if inequality of full dimensional cone is a
820 facet. If the routine returns true then the inequality is a
821 facet. If it returns false it is unknown.
822 */
fastIsFacetCriterion(IntegerVectorList const & normals,IntegerVectorList::const_iterator i)823 static bool fastIsFacetCriterion(IntegerVectorList const &normals, IntegerVectorList::const_iterator i)
824 {
825 int n=i->size();
826 for(int j=0;j<n;j++)
827 if((*i)[j]>0 || (*i)[j]<0)
828 {
829 int sign=(*i)[j]>0?1:-1;
830 bool isTheOnly=true;
831 for(IntegerVectorList::const_iterator k=normals.begin();k!=normals.end();k++)
832 if(k!=i)
833 {
834 if((*k)[j]*sign>0)
835 {
836 isTheOnly=false;
837 break;
838 }
839 }
840 if(isTheOnly)return true;
841 }
842 return false;
843 }
844
fastIsFacet(IntegerVectorList const & normals,IntegerVectorList::const_iterator i)845 bool fastIsFacet(IntegerVectorList const &normals, IntegerVectorList::const_iterator i)
846 {
847 if(fastIsFacetCriterion(normals,i))return true;
848 // log0 fprintf(Stderr,"LP\n");
849 return isFacet(normals,i);
850 }
851
852
853 /*
854 ONLY WORKS AFFINELY AND NOT PROJECTIVELY!!!
855 bool fastIsFacet(IntegerVectorList const &normals, IntegerVectorList::const_iterator i)
856 {
857 int n=i->size();
858
859 bool rightAnswer=isFacet(normals,i);
860
861 IntegerVectorList tempNormals=normals;
862 bool doLoop=true;
863 log0 AsciiPrinter(Stderr).printVector(*i);
864 do
865 {
866 log0 fprintf(Stderr,"TempNormal size:%i\n",tempNormals.size());
867 log0 AsciiPrinter(Stderr).printVectorList(tempNormals);
868 IntegerVector maxVector=*i;
869 IntegerVector minVector=*i;
870 for(IntegerVectorList::const_iterator k=tempNormals.begin();k!=tempNormals.end();k++)
871 {
872 maxVector=max(maxVector,*k);
873 minVector=min(minVector,*k);
874 }
875 IntegerVector maxAttained(n);
876 IntegerVector minAttained(n);
877 for(IntegerVectorList::const_iterator k=tempNormals.begin();k!=tempNormals.end();k++)
878 {
879 for(int j=0;j<n;j++)
880 {
881 if((*k)[j]==maxVector[j])maxAttained[j]++;
882 if((*k)[j]==minVector[j])minAttained[j]++;
883 }
884 }
885 int bestEntry=-1;
886 int bestCount=2000000000;
887 int bestValue=0;
888 for(int j=0;j<n;j++)
889 {
890 if((*i)[j]==maxVector[j])
891 {
892 if(maxAttained[j]<bestCount)
893 {
894 bestEntry=j;
895 bestValue=maxVector[j];
896 bestCount=maxAttained[j];
897 }
898 }
899 if((*i)[j]==minVector[j])
900 {
901 if(minAttained[j]<bestCount)
902 {
903 bestEntry=j;
904 bestValue=minVector[j];
905 bestCount=minAttained[j];
906 }
907 }
908 }
909 log0 fprintf(Stderr,"Best entry:%i, Best Count:%i bestValue : %i\n",bestEntry,bestCount,bestValue);
910
911 IntegerVectorList newList;
912 for(IntegerVectorList::const_iterator k=tempNormals.begin();k!=tempNormals.end();k++)
913 {
914 if((*k)[bestEntry]==bestValue)newList.push_back(*k);
915 }
916 if(newList.size()==1)
917 {
918 assert(rightAnswer==true);
919 return true;
920 }
921 doLoop=tempNormals.size()>newList.size();
922 tempNormals=newList;
923 }
924 while(doLoop);
925
926
927 log0 fprintf(Stderr,"LP\n");
928 return isFacet(normals,i);
929 }
930 */
931 /*
932 bool fastIsFacet(IntegerVectorList const &normals, IntegerVectorList::const_iterator i)
933 {
934 int n=i->size();
935
936 bool rightAnswer=isFacet(normals,i);
937
938 IntegerVectorList tempNormals=normals;
939 bool doLoop=true;
940 //log0 AsciiPrinter(Stderr).printVector(*i);
941 do
942 {
943 // log0 fprintf(Stderr,"TempNormal size:%i\n",tempNormals.size());
944 //log0 AsciiPrinter(Stderr).printVectorList(tempNormals);
945 IntegerVector maxVector=*i;
946 IntegerVector minVector=*i;
947 for(IntegerVectorList::const_iterator k=tempNormals.begin();k!=tempNormals.end();k++)
948 {
949 maxVector=max(maxVector,*k);
950 minVector=min(minVector,*k);
951 }
952
953
954
955 /* IntegerVector maxAttained(n);
956 IntegerVector minAttained(n);
957 for(IntegerVectorList::const_iterator k=tempNormals.begin();k!=tempNormals.end();k++)
958 {
959 for(int j=0;j<n;j++)
960 {
961 if((*k)[j]==maxVector[j])maxAttained[j]++;
962 if((*k)[j]==minVector[j])minAttained[j]++;
963 }
964 }
965 */
966 /* int bestEntry=-1;
967 int bestCount=2000000000;
968 int bestValue=0;
969 for(int j=0;j<n;j++)
970 {
971 if((*i)[j]==0)
972 {
973 if(maxVector[j]==0)
974 }
975 else
976 {
977 int sign=(*i)[j]>0?1:-1;
978
979 int numberOfVectorsWithSameSign=0;
980 for(IntegerVectorList::const_iterator k=normals.begin();k!=normals.end();k++)
981 {
982 if((*k)[j]*sign>0)
983 {
984 numberOfVectorsWithSameSign++;
985 }
986 }
987 if(numberOfVectorsWithSameSign==1)
988 return true;
989 }
990
991
992 if((*i)[j]==maxVector[j])
993 {
994 if(maxAttained[j]<bestCount)
995 {
996 bestEntry=j;
997 bestValue=maxVector[j];
998 bestCount=maxAttained[j];
999 }
1000 }
1001 if((*i)[j]==minVector[j])
1002 {
1003 if(minAttained[j]<bestCount)
1004 {
1005 bestEntry=j;
1006 bestValue=minVector[j];
1007 bestCount=minAttained[j];
1008 }
1009 }
1010 }
1011 //log0 fprintf(Stderr,"Best entry:%i, Best Count:%i bestValue : %i\n",bestEntry,bestCount,bestValue);
1012
1013 IntegerVectorList newList;
1014 for(IntegerVectorList::const_iterator k=tempNormals.begin();k!=tempNormals.end();k++)
1015 {
1016 if((*k)[bestEntry]==bestValue)newList.push_back(*k);
1017 }
1018 if(newList.size()==1)
1019 {
1020 assert(rightAnswer==true);
1021 return true;
1022 }
1023 doLoop=tempNormals.size()>newList.size();
1024 tempNormals=newList;
1025 }
1026 while(doLoop);
1027
1028
1029 // log0 fprintf(Stderr,"LP\n");
1030 return isFacet(normals,i);
1031 }
1032 */
1033
fastNormals(IntegerVectorList const & inequalities)1034 IntegerVectorList fastNormals(IntegerVectorList const &inequalities)
1035 {
1036 // log0 cerr << "Fast normals begin" << endl;
1037 // log0 fprintf(Stderr,"E");
1038 //log0 fprintf(Stderr,"Number of inequalities:%i\n",inequalities.size());
1039 IntegerVectorList normals=normalizedWithSumsAndDuplicatesRemoved(inequalities);
1040 // log0 fprintf(Stderr,"Number of inequalities:%i\n",normals.size());
1041 // log0 fprintf(Stderr,"F");
1042
1043
1044 // log0 AsciiPrinter(Stderr).printVectorList(normals);
1045 for(IntegerVectorList::iterator i=normals.begin();i!=normals.end();)
1046 //if(!isFacet(normals,i))
1047 if(!fastIsFacet(normals,i))
1048 {
1049 IntegerVectorList::iterator temp=i;
1050 i++;
1051 normals.erase(temp);
1052 }
1053 else i++;
1054
1055 // log0 fprintf(Stderr,"Number of inequalities:%i\n",normals.size());
1056
1057 // log0 fprintf(Stderr,"G");
1058 //gfan_log2 cerr << "Fast normals end" << endl;
1059 return normals;
1060 }
1061
1062 /* Virker ikke:
1063 IntegerVectorList fastNormals(IntegerVectorList const &inequalities)
1064 {
1065 log0 fprintf(Stderr,"E");
1066 log0 fprintf(Stderr,"Number of inequalities:%i\n",inequalities.size());
1067 IntegerVectorList normals=normalizedWithSumsAndDuplicatesRemoved(inequalities);
1068 log0 fprintf(Stderr,"Number of inequalities:%i\n",normals.size());
1069 log0 fprintf(Stderr,"F");
1070
1071
1072 // log0 AsciiPrinter(Stderr).printVectorList(normals);
1073
1074 bool didRemoveSomething;
1075 do
1076 {
1077 didRemoveSomething=false;
1078 for(IntegerVectorList::iterator i=normals.begin();i!=normals.end();i++)
1079 //if(!isFacet(normals,i))
1080 if(!fastIsFacetCriterion(normals,i))
1081 {
1082 IntegerVectorList::iterator temp=i;
1083 temp++;
1084 normals.erase(i);
1085 temp--;
1086 i=temp;
1087 didRemoveSomething=true;
1088 }
1089 }
1090 while(didRemoveSomething);
1091 for(IntegerVectorList::iterator i=normals.begin();i!=normals.end();i++)
1092 {
1093 if(!fastIsFacet(normals,i))
1094 {
1095 IntegerVectorList::iterator temp=i;
1096 temp++;
1097 normals.erase(i);
1098 temp--;
1099 i=temp;
1100 didRemoveSomething=true;
1101 }
1102 }
1103 log0 fprintf(Stderr,"Number of inequalities:%i\n",normals.size());
1104
1105 log0 fprintf(Stderr,"G");
1106 return normals;
1107 }
1108 */
1109
1110
normalConeOfMinkowskiSum(PolynomialSet const & polynomials,IntegerVector const & w)1111 PolyhedralCone normalConeOfMinkowskiSum(PolynomialSet const &polynomials, IntegerVector const &w)
1112 {
1113 // log0 cerr << "A";
1114 WeightReverseLexicographicTermOrder T(w);
1115 PolynomialSet g=polynomials;
1116 //log0 cerr << "I";
1117 g.markAndScale(T);
1118 //log0 cerr << "I";
1119 IntegerVectorList inequalities=fastNormals(wallInequalities(g));
1120 g=initialFormsAssumeMarked(g,w);
1121 IntegerVectorList equations=wallInequalities(g);
1122 //log0 cerr << "B";
1123 PolyhedralCone ret(inequalities,equations,g.getRing().getNumberOfVariables());
1124 //log0 cerr << "C";
1125 ret.canonicalize();
1126 //log0 cerr << "D";
1127 return ret;
1128 }
1129
1130
commonHomogeneitySpaceEquations(PolynomialSet const & polynomials)1131 IntegerVectorList commonHomogeneitySpaceEquations(PolynomialSet const &polynomials)
1132 {
1133 LexicographicTermOrder T;
1134 PolynomialSet g=polynomials;
1135 g.markAndScale(T);
1136 IntegerVectorList l=wallInequalities(g);
1137 FieldMatrix A=integerMatrixToFieldMatrix(rowsToIntegerMatrix(l,g.getRing().getNumberOfVariables()),Q);
1138 A.reduce();
1139 A.removeZeroRows();
1140 return fieldMatrixToIntegerMatrixPrimitive(A).getRows();
1141 }
1142
1143
commonHomogeneitySpaceGenerators(PolynomialSet const & polynomials)1144 IntegerVectorList commonHomogeneitySpaceGenerators(PolynomialSet const &polynomials)
1145 {
1146 LexicographicTermOrder T;
1147 PolynomialSet g=polynomials;
1148 g.markAndScale(T);
1149 IntegerVectorList l=wallInequalities(g);
1150 FieldMatrix A=integerMatrixToFieldMatrix(rowsToIntegerMatrix(l,g.getRing().getNumberOfVariables()),Q);
1151 A=A.reduceAndComputeKernel();
1152 return fieldMatrixToIntegerMatrixPrimitive(A).getRows();
1153 }
1154
1155
homogeneitySpaceEquations(Polynomial const & f)1156 FieldMatrix homogeneitySpaceEquations(Polynomial const &f)
1157 {
1158 assert(!f.isZero());
1159 IntegerVectorList l;
1160 TermMap::const_iterator j=f.terms.begin();
1161 IntegerVector first=j->first.exponent;
1162 for(j++;j!=f.terms.end();j++)
1163 l.push_back(first-j->first.exponent);
1164 return integerMatrixToFieldMatrix(rowsToIntegerMatrix(l,f.getRing().getNumberOfVariables()),Q);
1165 }
1166
homogeneitySpaceGenerators(Polynomial const & f)1167 FieldMatrix homogeneitySpaceGenerators(Polynomial const &f)
1168 {
1169 FieldMatrix A=homogeneitySpaceEquations(f);
1170 return A.reduceAndComputeKernel();
1171 }
1172
spanOfHomogeneitySpaces(PolynomialSet const & polynomials)1173 FieldMatrix spanOfHomogeneitySpaces(PolynomialSet const &polynomials)
1174 {
1175 FieldMatrix ret(Q,0,polynomials.getRing().getNumberOfVariables());
1176 for(PolynomialSet::const_iterator i=polynomials.begin();i!=polynomials.end();i++)
1177 ret=combineOnTop(ret,homogeneitySpaceGenerators(*i));
1178 ret.reduce();
1179 ret.removeZeroRows();
1180 return ret;
1181 }
1182