1 /*
2  * integergb.cpp
3  *
4  *  Created on: Dec 14, 2010
5  *      Author: anders
6  */
7 
8 #include "integergb.h"
9 #include "polynomial.h"
10 #include "field_rationals.h"
11 #include "printer.h"
12 #include "polyhedralcone.h"
13 #include "wallideal.h"
14 #include "tropical2.h"
15 
16 /*
17  * Implemented according to [Becker, Weispfenning] Chapter 10.1.
18  */
19 
dDivision(Polynomial p,PolynomialSet const & l,TermOrder const & termOrder)20 Polynomial dDivision(Polynomial p, PolynomialSet const &l, TermOrder const &termOrder)
21 {
22   PolynomialRing theRing=p.getRing();
23   Polynomial r(p.getRing());
24 
25   while(!p.isZero())
26     {
27       p.mark(termOrder);
28       Term initial=p.getMarked();
29       PolynomialSet::const_iterator i;
30       PolynomialSet::iterator j;
31       for(i=l.begin();i!=l.end();i++)
32         {
33           if(i->getMarked().m.exponent.divides(initial.m.exponent))
34             if((initial.c*(i->getMarked().c.inverse())).isInteger())break;
35         }
36 
37       {
38         if(i!=l.end())
39           {
40             Term s(-initial.c*i->getMarked().c.inverse(),Monomial(p.getRing(),initial.m.exponent-i->getMarked().m.exponent));
41             p.madd(s,*i);
42           }
43         else
44           {
45             p-=initial;
46             r+=initial;
47           }
48       }
49     }
50   return r;
51 }
52 
53 
eDivision(Polynomial p,PolynomialSet const & l,TermOrder const & termOrder)54 Polynomial eDivision(Polynomial p, PolynomialSet const &l, TermOrder const &termOrder)
55 {
56   PolynomialRing theRing=p.getRing();
57   Polynomial r(p.getRing());
58 
59 //  debug<<"eDivision input "<< p<<l<<termOrder;
60 
61   while(!p.isZero())
62     {
63       p.mark(termOrder);
64 //      debug<<"P "<<p<<"\n";
65       Term initial=p.getMarked();
66       PolynomialSet::const_iterator i;
67       PolynomialSet::iterator j;
68       for(i=l.begin();i!=l.end();i++)
69         {
70           if(i->getMarked().m.exponent.divides(initial.m.exponent))
71             {
72 //              debug<<"CHECKING:"<<initial<<"\n";
73               FieldElement q=integerDivision(initial.c,i->getMarked().c);
74               if(!q.isZero())break;
75             }
76         }
77 
78       {
79         if(i!=l.end())
80           {
81 //            debug<<"dividing by  "<<*i<<"\n";
82             Term s(-integerDivision(initial.c,i->getMarked().c),Monomial(p.getRing(),initial.m.exponent-i->getMarked().m.exponent));
83             p.madd(s,*i);
84           }
85         else
86           {
87             p-=initial;
88             r+=initial;
89           }
90       }
91     }
92   return r;
93 }
94 
sgpol(Polynomial const & g1,Polynomial const & g2,bool s)95 static Polynomial sgpol(Polynomial const &g1, Polynomial const &g2, bool s)
96 {
97   PolynomialRing R=g1.getRing();
98   FieldElement a1=g1.getMarked().c;
99   FieldElement a2=g2.getMarked().c;
100   IntegerVector t1=g1.getMarked().m.exponent;
101   IntegerVector t2=g2.getMarked().m.exponent;
102   FieldElement c1(Q);
103   FieldElement c2(Q);
104   FieldElement g=gcd(a1,a2,c1,c2);
105   FieldElement b1=a2*g.inverse();
106   FieldElement b2=a1*g.inverse();
107   IntegerVector s1=max(t1,t2)-t1;
108   IntegerVector s2=max(t1,t2)-t2;
109 
110 /*  debug
111   <<"a1 "<<a1<<"\n"
112   <<"a2 "<<a2<<"\n"
113 //  <<"t1 "<<t1<<"\n"
114 //  <<"t2 "<<t2<<"\n"
115   <<"c1 "<<c1<<"\n"
116   <<"c2 "<<c2<<"\n"
117   <<"g  "<<g<<"\n"
118   <<"b1 "<<b1<<"\n"
119   <<"b2 "<<b2<<"\n"
120   <<"s1 "<<s1<<"\n"
121   <<"s2 "<<s2<<"\n";
122 */
123   if(s)return Term(b1,Monomial(R,s1))*g1-Term(b2,Monomial(R,s2))*g2;
124   return Term(c1,Monomial(R,s1))*g1+Term(c2,Monomial(R,s2))*g2;
125 }
126 
127 
spol(Polynomial const & g1,Polynomial const & g2)128 Polynomial spol(Polynomial const &g1, Polynomial const &g2)
129 {
130   return sgpol(g1,g2,true);
131 }
132 
133 
gpol(Polynomial const & g1,Polynomial const & g2)134 Polynomial gpol(Polynomial const &g1, Polynomial const &g2)
135 {
136   return sgpol(g1,g2,false);
137 }
138 
zMinimize(PolynomialSet & F)139 void zMinimize(PolynomialSet &F)
140 {//can be done in nlogn time, but that requirese sorting according to coefficients.
141   for(PolynomialSet::iterator i=F.begin();i!=F.end();)
142     {
143       bool doDelete=false;
144       for(PolynomialSet::const_iterator j=F.begin();j!=F.end();j++)
145         if(i!=j)
146           if(j->getMarked().m.exponent.divides(i->getMarked().m.exponent))
147             if((j->getMarked().c.inverse()*i->getMarked().c).isInteger()){doDelete=true;break;}
148       if(doDelete)
149         {
150           PolynomialSet::iterator i2=i;
151           i++;
152           F.erase(i2);
153         }
154       else
155         i++;
156     }
157 }
158 
159 
zAutoReduce(PolynomialSet * g,TermOrder const & termOrder)160 void zAutoReduce(PolynomialSet *g, TermOrder const &termOrder)
161 {
162   for(PolynomialSet::iterator i=g->begin();i!=g->end();i++)
163     {
164 //      debug<<"1\n";
165       Polynomial temp(*i);
166       PolynomialSet::iterator tempIterator=i;
167       tempIterator++;
168       g->erase(i);
169       Monomial monomial=temp.getMarked().m;
170       if(temp.getMarked().c.sign()>0)
171         g->insert(tempIterator,eDivision(temp,*g,termOrder));
172       else
173         g->insert(tempIterator,eDivision(temp.getRing().zero()-temp,*g,termOrder));
174       tempIterator--;
175       i=tempIterator;
176       if(i->isZero())
177         {
178           assert(0);
179         }
180       else
181         i->mark(monomial);
182     }
183 }
184 
185 
186 
187 typedef pair<Polynomial,Polynomial> Pair;
zBuchberger(PolynomialSet & F,TermOrder const & T)188 void zBuchberger(PolynomialSet &F, TermOrder const &T)
189 {
190   list<Pair> B;
191   PolynomialSet G(F.getRing());
192   for(PolynomialSet::const_iterator i=F.begin();i!=F.end();i++)if(!i->isZero())G.push_back(*i);
193   G.mark(T);
194   for(PolynomialSet::const_iterator i=G.begin();i!=G.end();i++)
195     for(PolynomialSet::const_iterator j=G.begin();j!=i;j++)
196       B.push_back(Pair(*i,*j));
197   list<Pair> D;
198   list<Pair> C=B;
199 
200   while(!B.empty())
201     {
202 //      debug<<"Looping\n";
203 //      debug<<G;
204 
205       while(!C.empty())
206         {
207           Pair f=C.front();
208           C.pop_front();
209           PolynomialSet::const_iterator gi=G.begin();
210           f.first.mark(T);
211           f.second.mark(T);
212           IntegerVector lcm=max(f.first.getMarked().m.exponent,f.second.getMarked().m.exponent);
213           FieldElement HCf1=f.first.getMarked().c;
214           FieldElement HCf2=f.second.getMarked().c;
215           for(;gi!=G.end();gi++)
216             {
217               if(gi->getMarked().m.exponent.divides(lcm)
218                   && (gi->getMarked().c.inverse()*HCf1).isInteger()
219                   && (gi->getMarked().c.inverse()*HCf2).isInteger())break;
220             }
221           if(gi==G.end())
222             {
223   //            debug<<f.first.getMarked()<<"  "<<f.second.getMarked()<<"\n";
224 
225               Polynomial h=gpol(f.first,f.second);
226               Polynomial h0=dDivision(h,G,T);
227               if(h0.isZero())
228                 {
229   //                debug<<"remainder is zero!\n";
230                 }
231               else
232                 {
233                   h0.mark(T);
234                   for(PolynomialSet::const_iterator gi=G.begin();gi!=G.end();gi++)
235                     D.push_back(Pair(*gi,h0));
236                   G.push_back(h0);
237                 }
238             }
239         }
240       Pair f=B.front();
241       B.pop_front();
242       Polynomial h=spol(f.first,f.second);
243       Polynomial h0=dDivision(h,G,T);
244       h0.mark(T);
245       if(!h0.isZero())
246         {
247           for(PolynomialSet::const_iterator gi=G.begin();gi!=G.end();gi++)
248             D.push_back(Pair(*gi,h0));
249           G.push_back(h0);
250         }
251       C=D;
252       B.splice(B.end(),D);
253     }
254   F=G;
255   zMinimize(F);
256 //debug<<F<<"---------------------------------------------------------------------\n";
257   zAutoReduce(&F,T);
258 }
259 
260 
261 
262 
updatePolyhedralCone()263 void IntegerGroebnerFanTraverser::updatePolyhedralCone()
264 {
265   IntegerVectorList inequalities=fastNormals(wallInequalities(groebnerBasis));
266   IntegerVectorList empty;
267   theCone=PolyhedralCone(inequalities,empty,n);
268   theCone.canonicalize();
269  // debug<<theCone;
270  /// assert(0);
271 }
272 
IntegerGroebnerFanTraverser(PolynomialSet const & generators)273 IntegerGroebnerFanTraverser::IntegerGroebnerFanTraverser(PolynomialSet const &generators):
274 	ConeTraverser(generators.getRing().getNumberOfVariables()),
275   n(generators.getRing().getNumberOfVariables()),
276   groebnerBasis(generators)
277 {
278 //  LexicographicTermOrder tieBreaker;
279 //  zBuchberger(groebnerBasis,tieBreaker);
280 //debug<<"WE ASSUME THAT WE ALREADY HAVE A REDUCED GB\n";
281 
282   updatePolyhedralCone();
283 }
284 
changeCone(IntegerVector const & ridgeVector,IntegerVector const & rayVector)285 void IntegerGroebnerFanTraverser::changeCone(IntegerVector const &ridgeVector, IntegerVector const &rayVector)
286 {
287   IntegerVectorList tOld;
288   tOld.push_back(ridgeVector);
289   tOld.push_back(-rayVector);
290   MatrixTermOrder TOld(tOld);
291 
292   IntegerVectorList t;
293   t.push_back(ridgeVector);
294   t.push_back(rayVector);
295   MatrixTermOrder T(t);
296 
297 //debug<<"FLIP!\n";
298 //debug<<"OLDGB\n"<<groebnerBasis;
299 PolynomialSet g=initialFormsAssumeMarked(groebnerBasis,ridgeVector);
300 //debug<<"OLD INITIAL GB\n"<<g;
301 zBuchberger(g,T);
302 //debug<<"NEW INITIAL GB\n"<<g;
303 
304 
305 
306 PolynomialSet liftedBasis(groebnerBasis.getRing());
307 for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
308   {
309     Polynomial lt=*i;//Polynomial(i->getMarked());
310     Polynomial r=eDivision(lt, groebnerBasis,TOld);
311     Polynomial rIn=eDivision(lt,initialFormsAssumeMarked(groebnerBasis,ridgeVector),TOld);
312     liftedBasis.push_back(lt-r);
313 
314 //   debug<<"--------------------\n";
315 //    debug<<"lt:"<<lt<<"\n";
316 //    debug<<"r:"<<r<<"\n";
317 //    debug<<"rIn:"<<rIn<<"\n";
318 //    debug<<"--------------------\n";
319   }
320 liftedBasis.mark(T);
321 //debug<<"to be autoreduced:"<<liftedBasis;
322 //debug<<"autoreducing......\n";
323 liftedBasis.mark(T);
324 zAutoReduce(&liftedBasis,T);
325 PolynomialSet newInitialForms=initialForms(liftedBasis,ridgeVector);
326 //debug<<"LIFTED"<<liftedBasis;
327 //debug<<"NEWINITIALFORMS"<<newInitialForms;
328 
329 /*zBuchberger(groebnerBasis,T);
330 
331 {
332   groebnerBasis.sort_();
333   liftedBasis.sort_();
334   bool areEqual=(groebnerBasis.size()==liftedBasis.size());
335   if(areEqual)
336     {
337       PolynomialSet::const_iterator i=groebnerBasis.begin();
338       for(PolynomialSet::const_iterator j=liftedBasis.begin();j!=liftedBasis.end();j++,i++)
339         areEqual&=(*i-*j).isZero();
340     }
341   if(!areEqual)
342     {
343       debug<<groebnerBasis<<liftedBasis;
344       assert(0);
345     }
346 }*/
347 
348 groebnerBasis=liftedBasis;
349 groebnerBasis.mark(T);
350 
351 
352 //debug<<"NEW BASIS"<<groebnerBasis;
353 
354   updatePolyhedralCone();
355 }
356 
link(IntegerVector const & ridgeVector)357 IntegerVectorList IntegerGroebnerFanTraverser::link(IntegerVector const &ridgeVector)
358 {
359   IntegerVectorList ret;
360   IntegerVector v=theCone.link(ridgeVector).getUniquePoint();
361   ret.push_back(v);
362 
363   PolyhedralCone temp=intersection(PolyhedralCone::positiveOrthant(n),theCone.faceContaining(ridgeVector));
364   IntegerVector temp2=temp.getRelativeInteriorPoint();
365   if(temp2.min()>0)
366     {
367       ret.push_back(-v);
368     }
369 //  debug<<"LINK"<<ret;
370   return ret;
371 }
372 
373 
refToPolyhedralCone()374 PolyhedralCone & IntegerGroebnerFanTraverser::refToPolyhedralCone()
375 {
376   return theCone;
377 }
378