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