1 #include "buchberger.h"
2
3 #include <set>
4 #include <algorithm>
5 #include <iostream>
6
7 #include "division.h"
8 #include "printer.h"
9 #include "timer.h"
10 #include "parser.h"
11 #include "log.h"
12
13 #include "gebauermoeller.h"
14
15 static Timer buchbergerTimer("Buchberger",10);
16
17
sPolynomial(Polynomial a,Polynomial b)18 Polynomial sPolynomial(Polynomial a, Polynomial b)
19 {
20 bool comments=false;
21 if(comments)
22 {
23 AsciiPrinter(Stderr).printString("S(");
24 AsciiPrinter(Stderr).printPolynomial(a);
25 AsciiPrinter(Stderr).printString(",");
26 AsciiPrinter(Stderr).printPolynomial(b);
27 AsciiPrinter(Stderr).printString(")=");
28 }
29
30 //marked coefficient of a and b must be one
31
32 IntegerVector ina=a.getMarked().m.exponent;
33 IntegerVector inb=b.getMarked().m.exponent;
34
35 IntegerVector L=max(ina,inb);
36
37 if(comments)
38 AsciiPrinter(Stderr).printVector(L);
39
40
41 FieldElement const &f=a.getMarked().c;
42
43
44 Polynomial A=a;A*=Monomial(a.getRing(),L-ina);
45 Polynomial B=b;B*=Monomial(b.getRing(),L-inb);
46
47 if(comments)
48 {
49 AsciiPrinter(Stderr).printPolynomial(A-B);
50 AsciiPrinter(Stderr).printString("\n");
51 }
52
53 return A-B;
54 }
55
56
57 // Simple Buchberger
58
59 void buchbergerChain(PolynomialSet *g, TermOrder const &termOrder, bool allowSaturation);
60 void buchberger2(PolynomialSet *g, TermOrder const &termOrder, bool allowSaturation);
buchberger(PolynomialSet * g,TermOrder const & termOrder,bool allowSaturation)61 void buchberger/*Simple*/(PolynomialSet *g, TermOrder const &termOrder, bool allowSaturation)
62 {
63 int numberOfReductions;
64
65 return buchberger2(g, termOrder, allowSaturation);
66 //return buchbergerChain(g, termOrder, allowSaturation);
67 PolynomialRing theRing=g->getRing();
68 // gfan_log2 fprintf(Stderr,"ENTERING buchberger\n");
69 TimerScope ts(&buchbergerTimer);
70 PolynomialSet sPolynomials(theRing);
71
72 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
73 if(!i->isZero())sPolynomials.push_back(*i); // It is safe and useful to ignore the 0 polynomial
74
75 if(allowSaturation)sPolynomials.saturate();
76 sPolynomials.markAndScale(termOrder);
77
78 *g=PolynomialSet(theRing);
79
80 while(!sPolynomials.empty())
81 {
82 //pout<<int(sPolynomials.size())<<"\n";///////////////////////
83 Polynomial p=*sPolynomials.begin();
84 sPolynomials.pop_front();
85
86 p=division(p,*g,termOrder);
87 numberOfReductions++;
88 if(!p.isZero())
89 {
90 // pout<<"NONZERO\n";
91 if(allowSaturation)p.saturate();
92 p.mark(termOrder);
93 p.scaleMarkedCoefficientToOne();
94 bool isMonomial=p.isMonomial();
95 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
96 if((!isMonomial) || (!i->isMonomial())) // 2 % speed up!
97 {
98 if(!relativelyPrime(i->getMarked().m.exponent,p.getMarked().m.exponent))
99 {
100 Polynomial s=sPolynomial(*i,p);
101 s.mark(termOrder); // with respect to some termorder
102 s.scaleMarkedCoefficientToOne();
103 sPolynomials.push_back(s);
104 }
105 }
106
107 if(1){//Mora trick (from his book). This is not exactly what Mora suggests, but it seems to work. (Book1, Part 3, Remark 3.6.1)
108 for(PolynomialSet::iterator i=g->begin();i!=g->end();)
109 if(p.getMarked().m.exponent.divides(i->getMarked().m.exponent))
110 {
111 PolynomialSet::iterator j=i;j++;
112 Polynomial q=sPolynomial(*i,p);
113 if(allowSaturation)q.saturate();
114 if(!q.isZero())
115 {
116 sPolynomials.push_back(q);
117 }
118 g->erase(i);
119 // debug<<"erasing\n";
120
121 i=j;
122 }
123 else
124 i++;
125 }
126
127
128
129 g->push_back(p);
130 {
131 static int t;
132 t++;
133 // if((t&31)==0)fprintf(Stderr," gsize %i spolys:%i\n",g->size(),sPolynomials.size());
134 }
135 }
136 // else
137 // pout<<"ZERO\n";
138 }
139 //gfan_log2 fprintf(Stderr," buchberger minimize\n");
140 minimize(g);
141 //gfan_log2 fprintf(Stderr," buchberger autoreduce\n");
142 autoReduce(g,termOrder);
143 //gfan_log2 fprintf(Stderr,"LEAVING buchberger\n\n");
144
145 cerr<<"NumberOfReductions: "<<numberOfReductions<<std::endl;
146 }
147
148
149
150 // Buchberger with chain criterion
151
152 struct ChainPair
153 {
154 PolynomialSet::const_iterator a,b;
155 int A,B;
156 IntegerVector lcm;
ChainPairChainPair157 ChainPair(PolynomialSet::const_iterator const &a_,PolynomialSet::const_iterator const &b_,int A_,int B_):
158 a(a_),
159 b(b_),
160 A(A_),
161 B(B_),
162 lcm(max(a_->getMarked().m.exponent,b_->getMarked().m.exponent))
163 {
164 if(B<A)
165 {
166 a=b_;
167 b=a_;
168 A=B_;
169 B=A_;
170 }
171 }
operator <ChainPair172 bool operator<(const ChainPair & b)const
173 {
174 if(b.lcm.sum()<lcm.sum())return false;
175 if(lcm.sum()<b.lcm.sum())return true;
176 if(b.lcm<lcm)return true;
177 if(lcm<b.lcm)return false;
178 if(A<b.A)return true;
179 if(A>b.A)return false;
180 if(B<b.B)return true;
181 if(B>b.B)return false;
182 return false;
183 // assert(0);
184 }
185 };
186
187 typedef set<ChainPair> ChainPairList;
188
canBeRemoved(ChainPairList const & P,IntegerVector const & lcm,PolynomialSet::const_iterator i,PolynomialSet::const_iterator l,int I,int L)189 static bool canBeRemoved(ChainPairList const &P, IntegerVector const &lcm, PolynomialSet::const_iterator i, PolynomialSet::const_iterator l, int I, int L)
190 {
191 for(ChainPairList::const_iterator t=P.begin();t!=P.end();t++)
192 {
193 if(t->a==i && t->b!=l && t->b->getMarked().m.exponent.divides(lcm))
194 {
195 if(P.count(ChainPair(l,t->b,L,t->B))==1)return true;
196 }
197 if(t->b==i && t->a!=l && t->a->getMarked().m.exponent.divides(lcm))
198 {
199 if(P.count(ChainPair(l,t->a,L,t->A))==1)return true;
200 }
201 if(t->a==l && t->b!=i && t->b->getMarked().m.exponent.divides(lcm))
202 {
203 if(P.count(ChainPair(i,t->b,I,t->B))==1)return true;
204 }
205 if(t->b==l && t->a!=i && t->a->getMarked().m.exponent.divides(lcm))
206 {
207 if(P.count(ChainPair(i,t->a,I,t->A))==1)return true;
208 }
209 }
210 return false;
211
212 // return false;
213 for(ChainPairList::const_iterator t=P.begin();t!=P.end();t++)
214 {
215 if(t->a==i && t->b!=l && t->b->getMarked().m.exponent.divides(lcm) /*||
216 t->b==i && t->a!=l && t->a->getMarked().m.exponent.divides(lcm) ||
217 t->a==l && t->b!=i && t->b->getMarked().m.exponent.divides(lcm) ||
218 t->b==l && t->a!=i && t->a->getMarked().m.exponent.divides(lcm)*/)return true;
219 }
220 return false;
221 }
222
printPairs(ChainPairList const & P)223 void printPairs(ChainPairList const &P)
224 {
225 // return;
226 for(ChainPairList::const_iterator t=P.begin();t!=P.end();t++)
227 {
228 cerr<<"("<<t->A<<","<<t->B<<")[";
229 AsciiPrinter(Stderr)<<t->a->getMarked()<<","<<t->b->getMarked()<<"]<"<<t->lcm <<">";
230 cerr<<endl;
231 }
232 cerr<<endl;
233 }
234
235 //void buchbergerChain(PolynomialSet *g, TermOrder const &termOrder)
buchbergerChain(PolynomialSet * g,TermOrder const & termOrder,bool allowSaturation)236 void buchbergerChain(PolynomialSet *g, TermOrder const &termOrder, bool allowSaturation)
237 {
238 int nDivisions=0;
239 int nZeroDivisions=0;
240 int nRemovedChain=0;
241 PolynomialRing theRing=g->getRing();
242 TimerScope ts(&buchbergerTimer);
243 //cerr<<g->size();
244 {
245 PolynomialSet g2(theRing);
246 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
247 if(!i->isZero())g2.push_back(*i); // It is safe and useful to ignore the 0 polynomial
248
249 *g=g2;
250 if(allowSaturation)g->saturate();
251 }
252 g->markAndScale(termOrder);
253
254 ChainPairList P;//use better data structure for this
255 ChainPairList Pchecked;//use better data structure for this
256 int I=0;
257 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++,I++)
258 {
259 int J=0;
260 for(PolynomialSet::const_iterator j=g->begin();j!=i;j++,J++)
261 {
262 if(relativelyPrime(i->getMarked().m.exponent,j->getMarked().m.exponent) || (i->isMonomial() && j->isMonomial()))
263 Pchecked.insert(ChainPair(j,i,J,I));//here
264 else
265 P.insert(ChainPair(j,i,J,I));//here
266 }
267 }
268
269 while(!P.empty())
270 {
271 // pout<<*g;
272 // cerr<<"P\n";printPairs(P);cerr<<"Pchecked\n";printPairs(Pchecked);
273
274
275 PolynomialSet::const_iterator i=P.begin()->a;
276 PolynomialSet::const_iterator l=P.begin()->b;
277 int I=P.begin()->A;
278 int L=P.begin()->B;
279
280 if(relativelyPrime(i->getMarked().m.exponent,l->getMarked().m.exponent) || (i->isMonomial() && l->isMonomial()))
281 {
282 // Pchecked.push_back(P.front());
283 // P.pop_front();
284 Pchecked.insert(*P.begin());
285 P.erase(P.begin());
286 }
287 else
288 {
289 IntegerVector lcm=max(i->getMarked().m.exponent,l->getMarked().m.exponent);
290
291 if(/*canBeRemoved(P,lcm,i,l) ||*/ canBeRemoved(Pchecked,lcm,i,l,I,L))
292 {
293 //cerr<<"removin"<<endl;
294 //P.pop_front(); // This might remove elements from P with a trivial S-poly
295 P.erase(P.begin()); // This might remove elements from P with a trivial S-poly
296 nRemovedChain++;
297 }
298 else
299 {
300 Polynomial p=sPolynomial(*i,*l);
301
302 p.mark(termOrder);
303 p=division(p,*g,termOrder);nDivisions++;
304 if(allowSaturation)p.saturate();
305 // Pchecked.push_back(P.front());
306 // P.pop_front();
307 Pchecked.insert(*P.begin());
308 P.erase(P.begin());
309
310 if(!p.isZero())
311 {
312 p.mark(termOrder);
313 p.scaleMarkedCoefficientToOne();
314 int K=g->size();
315 g->push_back(p);
316 PolynomialSet::const_iterator k=g->end();k--;
317 int I=0;
318 for(PolynomialSet::const_iterator i=g->begin();i!=k;i++,I++)
319 {
320 if(relativelyPrime(i->getMarked().m.exponent,k->getMarked().m.exponent) || (i->isMonomial() && k->isMonomial()))
321 Pchecked.insert(ChainPair(i,k,I,K));//here
322 else
323 P.insert(ChainPair(i,k,I,K));//here
324 // P.insert(ChainPair(i,k,I,K));//here
325 }
326 }
327 else
328 nZeroDivisions++;
329 }
330 }
331 }
332 // AsciiPrinter(Stderr)<<*g;
333 minimize(g);
334 autoReduce(g,termOrder);
335 cerr<<"NDIV "<<nDivisions<<" NREMOVEDCHAIN "<<nRemovedChain<<" NZERODIVISIONS "<<nZeroDivisions<<endl;
336 }
337
338
339
340 /*class MyPolynomialCompare
341 {
342 TermOrder const &termOrder;
343 public:
344 MyPolynomialCompare(TermOrder const &termOrder_):termOrder(termOrder_)
345 {
346
347 }
348 bool operator()(const Polynomial &a, const Polynomial &b)const
349 {
350 if(termOrder(a.getMarked().m.exponent,b.getMarked().m.exponent)<0)return true;
351 if(termOrder(a.getMarked().m.exponent,b.getMarked().m.exponent)>0)return false;
352 return PolynomialCompare(a,b);
353 }
354 };*/
355
356
autoReduceNonGB(PolynomialSet * g,TermOrder const & termOrder)357 void autoReduceNonGB(PolynomialSet *g, TermOrder const &termOrder)
358 {
359 g->markAndScale(termOrder);
360 // g->sort(PolynomialCompareMarkedTerms(termOrder));
361 g->sort(PolynomialCompareMarkedTerms(TotalDegreeTieBrokenTermOrder(termOrder)));
362 // debug<<"Sorted:"<<*g;
363
364 for(PolynomialSet::iterator i=g->begin();i!=g->end();)
365 {
366 Polynomial temp(*i);
367 PolynomialSet::iterator tempIterator=i;
368 tempIterator++;
369 g->erase(i);
370 Monomial monomial=temp.getMarked().m;
371
372 Polynomial remainder=division(temp,*g,termOrder);
373 // debug<<remainder<<"\n";
374 if(!remainder.isZero())
375 {
376 g->insert(tempIterator,remainder);
377 //g->insert(tempIterator,smartDivision(temp,*g,termOrder));
378 tempIterator--;
379 i=tempIterator;
380 i->mark(termOrder);
381 // i->mark(monomial);
382 i++;
383 }
384 else
385 i=tempIterator;
386 }
387 }
388
389 struct polynomialOrder
390 {
391 TermOrder const ℴ
polynomialOrderpolynomialOrder392 polynomialOrder(TermOrder const &order_):order(order_){}
393
operator ()polynomialOrder394 bool operator()(const Polynomial &a, const Polynomial &b)
395 {
396 return order(a.getMarked().m.exponent, b.getMarked().m.exponent);
397 }
398 };
399
400 /*
401 * g must be marked.
402 * all initial monomials in g must be different
403 */
autoReduceSorting(PolynomialSet * g,TermOrder const & termOrder)404 void autoReduceSorting(PolynomialSet *g, TermOrder const &termOrder)
405 {
406 g->sort(polynomialOrder(termOrder));
407 PolynomialSet temp(g->getRing());
408 for(PolynomialSet::iterator i=g->begin();i!=g->end();i++)
409 {
410 Monomial monomial=i->getMarked().m;
411 Polynomial f=division(*i,temp,termOrder);
412 f.mark(monomial);
413 temp.push_back(f);
414 }
415 *g=temp;
416 }
417
418 /* Mora's example from his book:
419 Q[v,w,z,y,x]
420 {v2-xz,y2-x3,yzv-x2w}
421 */
buchberger2(PolynomialSet * g,TermOrder const & termOrder,bool allowSaturation)422 void buchberger2(PolynomialSet *g, TermOrder const &termOrder, bool allowSaturation)
423 {
424 bool printDebug=false;
425 PolynomialRing theRing=g->getRing();
426
427 // debug<<"Buchberger2 on:"<<theRing<<*g<<"\n";
428
429
430 g->markAndScale(termOrder);
431 autoReduceNonGB(g,termOrder);
432 g->markAndScale(termOrder);
433 g->sort(PolynomialCompareMarkedTermsReverse(termOrder));
434 // g->computeInitialSugar();
435
436 if(allowSaturation)
437 {
438 for(PolynomialSet::iterator i=g->begin();i!=g->end();i++)
439 i->saturate();
440 }
441
442 // debug<<"Reduced NonGB:"<<theRing<<*g<<"\n";
443
444 // SPairContainerPair sPolynomials(theRing,termOrder);
445 SPairContainerOptimized sPolynomials(theRing,termOrder);
446 // SPairContainerOptimizedPacked sPolynomials(theRing,termOrder,g);
447
448 vector<Polynomial> G;
449 int indexj=0;
450 for(PolynomialSet::iterator j=g->begin();j!=g->end();j++,indexj++)
451 {
452 G.push_back(*j);
453 sPolynomials.updatePairs(G,indexj);
454 if(printDebug)
455 {
456 // printSPairSet(sPolynomials);
457 cout<<"-----\n";
458 }
459 }
460
461 // printSPairSet(sPolynomials);
462
463 int numberOfCriticalPairsConsidered=0;
464 int numberOfUsefulCriticalPairs=0;
465
466 while(!sPolynomials.isEmpty())
467 {
468 pair<int,int> ij=sPolynomials.popPair();
469 Polynomial p=sPolynomial(G[ij.first],G[ij.second]);
470 // if(printDebug)
471 // cout<<"Processing"<<sPolynomials.begin()->i+1<<sPolynomials.begin()->j+1<<"\n";
472 p.mark(termOrder);
473 p.scaleMarkedCoefficientToOne();
474
475 numberOfCriticalPairsConsidered++;
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490 p=division(p,*g,termOrder);
491 if(allowSaturation)p.saturate();
492 if(!p.isZero())
493 {
494 p.mark(termOrder);
495 p.scaleMarkedCoefficientToOne();
496 g->push_back(p);
497 G.push_back(p);
498 numberOfUsefulCriticalPairs++;
499 gfan_log2
500 {
501 static int t;
502 if(((++t)&=31)==0)
503 fprintf(Stderr,"gsize:%i spolys:%i n:%i\n",(int)g->size()+1,(int)sPolynomials.size(),(int)g->getRing().getNumberOfVariables());
504 }
505 sPolynomials.updatePairs(G,indexj);//,truncationDegree,&grading);
506
507 indexj++;
508 if(printDebug)
509 {
510 //printSPairSet(sPolynomials);
511 cout<<"-----\n";
512 }
513 }
514 else
515 if(printDebug) cout<<"zero reduction\n";
516 }
517 *g=PolynomialSet(theRing);
518 for(PolynomialVector::const_iterator i=G.begin();i!=G.end();i++)g->push_back(*i);
519 // debug<<"MINIMIZING\n";
520 minimize(g);
521 // debug<<"AUTOREDUCING\n";
522 // debug.printPolynomialSet(*g,true);
523 autoReduceSorting(g,termOrder);
524 // debug.printPolynomialSet(*g,true);
525 // autoReduce(g,termOrder);
526 // debug<<"RETURNING\n";
527
528 // if(printComments)
529 // fprintf(Stderr,"Number of critical pairs considered(divisions) %i Useful %i\n",numberOfCriticalPairsConsidered,numberOfUsefulCriticalPairs);
530 }
531
532
533
minimize(PolynomialSet * g)534 void minimize(PolynomialSet *g)
535 {//CHECK THAT THIS ROUTINE WORKS IF TWO GENERATORS HAVE THE SAME INITIAL TERM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
536 PolynomialRing theRing=g->getRing();
537 PolynomialSet ret(theRing);
538
539 g->sort(polynomialOrder(LexicographicTermOrder()));
540
541 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
542 {
543 bool someDivides=false;
544 for(PolynomialSet::const_iterator j=ret.begin();j!=ret.end();j++)
545 // for(PolynomialSet::const_iterator j=g->begin();j!=i;j++) //changed Feb 2009
546 {
547 if(j->getMarked().m.exponent.divides(i->getMarked().m.exponent))
548 {
549 someDivides=true;
550 }
551 }
552 if(!someDivides)
553 ret.push_back(*i);
554 }
555
556 *g=ret;
557 }
558
559
560
autoReduce(PolynomialSet * g,TermOrder const & termOrder)561 void autoReduce(PolynomialSet *g, TermOrder const &termOrder)
562 {
563 /**
564 * TODO: there should be two options : supplying a termorder, or not supplying a termorder. In the latter case this routine should decide if it wants to compute one.
565 */
566 // WeightTermOrder termOrder2(termorderWeight(*g));//REMOVE ME ?? JAN 2009
567
568
569 for(PolynomialSet::iterator i=g->begin();i!=g->end();i++)
570 {
571 Polynomial temp(*i);
572 PolynomialSet::iterator tempIterator=i;
573 tempIterator++;
574 g->erase(i);
575 Monomial monomial=temp.getMarked().m;
576 g->insert(tempIterator,division(temp,*g,termOrder));
577 //g->insert(tempIterator,smartDivision(temp,*g,termOrder));
578 tempIterator--;
579 i=tempIterator;
580 i->mark(monomial);
581 }
582 }
583
584
isMarkedGroebnerBasis(PolynomialSet const & g)585 bool isMarkedGroebnerBasis(PolynomialSet const &g)
586 {
587 int counter=0;
588 for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
589 {
590 gfan_log2 fprintf(Stderr,"%i ",counter++);
591 for(PolynomialSet::const_iterator j=i;j!=g.end();j++)
592 if(!relativelyPrime(i->getMarked().m.exponent,j->getMarked().m.exponent))
593 {
594 Polynomial s=sPolynomial(*i,*j);
595 if(!division(s,g,LexicographicTermOrder()).isZero())
596 {
597 log3{AsciiPrinter(Stderr)<<"Spoly("<<*i<<","<<*j<<")="<<sPolynomial(*i,*j)<<" gives remainder "<< division(s,g,LexicographicTermOrder()) <<"\n";}
598 return false;
599 }
600 }
601 }
602 return true;
603 }
604
605
isMinimal(PolynomialSet const & g)606 bool isMinimal(PolynomialSet const &g)
607 {
608 PolynomialSet temp=g.markedTermIdeal();
609 minimize(&temp);
610 return temp.size()==g.size();
611 }
612
613
isReduced(PolynomialSet const & g)614 bool isReduced(PolynomialSet const &g)
615 {
616 if(!isMinimal(g))return false;
617 PolynomialSet temp=g.markedTermIdeal();
618 for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
619 for(TermMap::const_iterator j=i->terms.begin();j!=i->terms.end();j++)
620 if(!(j->first.exponent==i->getMarked().m.exponent))
621 for(PolynomialSet::const_iterator k=temp.begin();k!=temp.end();k++)
622 if(k->getMarked().m.exponent.divides(j->first.exponent))return false;
623 return true;
624 }
625