1 #include "tropicalbasis.h"
2
3 #include <iostream>
4
5 #include "buchberger.h"
6 #include "groebnerengine.h"
7 #include "tropical.h"
8 #include "tropical2.h"
9 #include "division.h"
10 #include "wallideal.h"
11 #include "halfopencone.h"
12 #include "log.h"
13
14 #include "timer.h"
15
16 static Timer iterativeTropicalBasisTimer("Iterative tropical basis",1);
17
18 typedef set<IntegerVector> IntegerVectorSet;
19
cleverSaturation(Polynomial const & p,IntegerVector const & w)20 static Polynomial cleverSaturation(Polynomial const &p, IntegerVector const &w)
21 {
22 PolynomialRing theRing=p.getRing();
23 if(p.isZero())return p;
24 if(w.size()==0)return p;
25 Polynomial f=initialForm(p,w);
26 debug<<"P:"<<p<<" w: "<<w<<" f: "<<f<<"\n";
27 IntegerVector gcd=f.greatestCommonMonomialDivisor();
28 if(!gcd.isZero())
29 {
30 debug<<"OLD"<<p<<"\n";
31 debug<<"initialForm"<<initialForm(p,w)<<"\n";
32 PolynomialSet reducer(p.getRing());
33 reducer.push_back(theRing.monomialFromExponentVector(IntegerVector::allOnes(theRing.getNumberOfVariables()))-theRing.one());
34 reducer.markAndScale(LexicographicTermOrder());
35 debug<<gcd;
36 int s=gcd.max();
37 Polynomial all=theRing.monomialFromExponentVector(s*IntegerVector::allOnes(theRing.getNumberOfVariables())-gcd);
38 Polynomial g= division(all*p,reducer,LexicographicTermOrder());
39 g.saturate();
40 debug<<"NEW"<<g<<"\n";
41 return g;
42 }
43 return p;
44 }
45
initialSaturatingBuchberger(PolynomialSet * g,TermOrder const & termOrder_,IntegerVector const & w_)46 static void initialSaturatingBuchberger(PolynomialSet *g, TermOrder const &termOrder_, IntegerVector const &w_)
47 {
48 int n=w_.size();
49 IntegerVectorList temp;
50 temp.push_back(IntegerVector::standardVector(n+1,n));
51 IntegerVector w=concatenation(w_,IntegerVector(1));w[w.size()-1]=-w_.sum();
52 temp.push_back(w);
53 MatrixTermOrder termOrder(temp);
54 PolynomialRing theRing=g->getRing().withVariablesAppended("Z");
55 PolynomialSet sPolynomials(theRing);
56
57 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
58 if(!i->isZero())sPolynomials.push_back(cleverSaturation(i->embeddedInto(theRing),w)); // It is safe and useful to ignore the 0 polynomial
59
60 sPolynomials.push_back(theRing.monomialFromExponentVector(IntegerVector::allOnes(theRing.getNumberOfVariables()))-theRing.one());
61
62 sPolynomials.saturate();
63 sPolynomials.markAndScale(termOrder);
64
65 *g=PolynomialSet(theRing);
66
67 while(!sPolynomials.empty())
68 {
69 Polynomial p=*sPolynomials.begin();
70 sPolynomials.pop_front();
71
72 p=division(p,*g,termOrder);
73 if(!p.isZero())
74 {
75 p=cleverSaturation(p,w);
76 p.mark(termOrder);
77 p.scaleMarkedCoefficientToOne();
78 bool isMonomial=p.isMonomial();
79 for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
80 if((!isMonomial) || (!i->isMonomial())) // 2 % speed up!
81 {
82 if(!relativelyPrime(i->getMarked().m.exponent,p.getMarked().m.exponent))
83 {
84 Polynomial s=sPolynomial(*i,p);
85 s.mark(termOrder); // with respect to some termorder
86 s.scaleMarkedCoefficientToOne();
87 sPolynomials.push_back(s);
88 }
89 }
90 g->push_back(p);
91 {
92 static int t;
93 t++;
94 // if((t&31)==0)fprintf(Stderr," gsize %i spolys:%i\n",g->size(),sPolynomials.size());
95 }
96 }
97 }
98 minimize(g);
99 autoReduce(g,termOrder);
100 }
101
102
103
tropicalBasisOfCurve(int n,PolynomialSet g,PolyhedralFan * intersectionFan,int linealitySpaceDimension)104 PolynomialSet tropicalBasisOfCurve(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension) //Assuming g is homogeneous
105 {
106
107 /*
108 * TODO: Introduce the notion of a tropical prebasis:
109 *
110 * Definition. A set of polynomials f_1,...,f_m is called a tropical prebasis for the ideal they
111 * generate if for every w not in the tropical variety of that ideal there exists a monomial in
112 * the ideal generated by the initial forms of f_1,...,f_m w.r.t. w.
113 *
114 * Computing a tropical prebasis could be faster than computing a tropical basis since fewer
115 * groebner bases for the originial ideal might be needed. Still, however, it is relatively easy
116 * to determine the tropical variety given a tropical prebasis.
117 */
118 PolynomialSet originalG=g;
119 // bool prebasis=true;
120 // debug<<"PREBASIS="<<prebasis<<"\n";
121 log2 debug<<"TropicalBasis begin\n";
122 log2 debug<<g;
123 int homog=linealitySpaceDimension;
124 log2 D(linealitySpaceDimension);
125 assert(homog>0 || n==0);
126 TimerScope ts(&iterativeTropicalBasisTimer);
127 PolyhedralFan f(n);
128 if(!intersectionFan)intersectionFan=&f;
129
130 // *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);
131 // log1 fprintf(Stderr,"WARINING USING EXPERIMENTAL TROPICAL HYPERSURFACE INTERSECTION ROUTINE!!\n");
132 *intersectionFan=tropicalHyperSurfaceIntersectionClosed(n, g);
133
134 IntegerVectorSet containsNoMonomialCache;
135
136 while(1)
137 {
138 PolyhedralFan::coneIterator i;
139
140 // {AsciiPrinter P(Stderr);intersectionFan->printWithIndices(&P);}
141 restart:
142 // {AsciiPrinter P(Stderr);intersectionFan->printWithIndices(&P);}
143 for(i=intersectionFan->conesBegin();i!=intersectionFan->conesEnd();i++)
144 {
145 // log1 cerr<<"!@#$";
146 IntegerVector relativeInteriorPoint=i->getRelativeInteriorPoint();
147 // log1 cerr<<"1234/n";
148
149 if(containsNoMonomialCache.count(relativeInteriorPoint)>0)
150 {
151 log2 fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
152 }
153 else
154 {
155 /* if(prebasis)
156 {
157 if(containsMonomial(initialForms(g,relativeInteriorPoint)))
158 {
159 intersectionFan->insertFacetsOfCone(*i);
160 intersectionFan->remove(*i);
161 debug<<"LOWERING DIMENSION OF CONE\n";//TODO: checking cones in order of dimension could avoid this.
162 goto restart;
163 }
164 }*/
165 WeightReverseLexicographicTermOrder t(relativeInteriorPoint);
166 log2 fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
167 log2 AsciiPrinter(Stderr).printVector(relativeInteriorPoint);
168 log2 fprintf(Stderr,"\n");
169 PolynomialSet h2=originalG;//g;//<------------------- Using the original set here speeds up things a lot in the starting cone via stable intersections algorithm
170 // debug<<"g"<<g;
171
172 /* {Adjust these lines somehow to enable the saturating buchberger.
173 debug<<h2;
174 initialSaturatingBuchberger(&h2, t, relativeInteriorPoint);
175 //buchberger(&h2,t,true);
176 debug<<"SATURATING BUCHBERGER DONE\n";
177 }*/
178 // buchberger(&h2,t);
179 h2=GE_groebnerBasis(h2,t,true/*autoreduce*/,true/*saturate*/);
180 // debug<<"h2"<<h2;
181 log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
182
183 // debug<<h2;
184 // log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
185
186 PolynomialSet wall=initialFormsAssumeMarked(h2,relativeInteriorPoint);
187
188 if(containsMonomial(wall))
189 {
190 log2 fprintf(Stderr,"Initial ideal contains a monomial.\n");
191 Polynomial m(computeTermInIdeal(wall));
192 log2 fprintf(Stderr,"Done computing term in ideal\n");
193
194 // Polynomial temp=m-division(m,h2,LexicographicTermOrder());
195 Polynomial temp=m-division(m,h2,t);
196 g.push_back(temp);
197
198 log2 fprintf(Stderr,"Adding element to basis:\n");
199 log2 AsciiPrinter(Stderr).printPolynomial(temp);
200 log2 fprintf(Stderr,"\n");
201
202 *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
203 break;
204 }
205 else
206 {
207 if(i->dimension()<=1+homog)
208 //if(!containsMonomial(wall) && i->dimension()<=1+homog)//line for testing perturbation code
209 {
210 log2 fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
211 containsNoMonomialCache.insert(relativeInteriorPoint);
212 }
213 else
214 {
215 /* We need to compute the initial ideal with
216 respect to "relativeInteriorPoint" perturbed
217 with a basis of the span of the cone. Instead
218 of perturbing we may as well compute initial
219 ideal successively. We have already computed
220 the initial ideal with respect to
221 "relativeInteriorPoint". To get the perturbed
222 initial ideal we take initial ideal with
223 respect to each vector in the basis of the
224 span.*/
225 IntegerVectorList empty;
226 PolyhedralCone dual=PolyhedralCone(empty,i->getEquations(),i->ambientDimension()).dualCone();
227 dual.canonicalize();
228 IntegerVectorList basis=dual.getEquations();
229 PolynomialSet witnessLiftBasis=h2;//basis with respect to relativeInteriorPoint
230 log2 debug<<"basis"<<basis<<"\n";
231 for(IntegerVectorList::const_iterator j=basis.begin();j!=basis.end();j++)
232 {
233 log2 debug<<"wall"<<wall<<"\n";
234 WeightReverseLexicographicTermOrder t(*j);
235 PolynomialSet h3=wall;
236 // buchberger(&h3,t);
237 h3=GE_groebnerBasis(h3,t,true/*autoreduce*/,false/*saturate*/);
238 wall=initialFormsAssumeMarked(h3,*j);
239 witnessLiftBasis=liftBasis(h3,witnessLiftBasis);
240 }
241 log2 debug<<"wall"<<wall<<"\n";
242 if(containsMonomial(wall))
243 {
244 Polynomial m(computeTermInIdeal(wall));
245 Polynomial temp=m-division(m,witnessLiftBasis,LexicographicTermOrder());
246 g.push_back(temp);
247 *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
248 break;
249 }
250 else
251 {
252 debug<<"ERROR: INITIAL IDEAL DOES NOT CONTAIN A MONOMIAL AS EXPECTED.\nLIKELY CAUSE: THE USER DID NOT SPECIFY AN IDEAL SATISFYING THE REQUIREMENT THAT IT SHOULD DEFINE A \"TROPICAL CURVE\".";
253 assert(0);
254 }
255 }
256 }
257 }
258 }
259 if(i==intersectionFan->conesEnd())break;
260 }
261
262 log2 debug<<"TropicalBasis end\n";
263 log2 cerr <<"RETURNING";
264 return g;
265 }
266
267 /*
268 PolynomialSet iterativeTropicalBasisNoPerturbation(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension, bool doPrint) //Assuming g is homogeneous
269 {
270 TimerScope ts(&iterativeTropicalBasisTimer);
271 PolyhedralFan f(n);
272 if(!intersectionFan)intersectionFan=&f;
273
274 *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);
275
276 IntegerVectorSet containsNoMonomialCache;
277
278 while(1)
279 {
280 // AsciiPrinter(Stderr).printPolyhedralFan(*intersectionFan);
281 // assert(f.getMaxDimension()==1);
282
283 IntegerVectorList l=intersectionFan->getRelativeInteriorPoints();
284
285 IntegerVectorList::const_iterator i;
286 for(i=l.begin();i!=l.end();i++)
287 {
288 if(containsNoMonomialCache.count(*i)>0)
289 {
290 if(doPrint)fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
291 }
292 else
293 {
294 WeightReverseLexicographicTermOrder t(*i);
295 if(doPrint)fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
296 if(doPrint)AsciiPrinter(Stderr).printVector(*i);
297 if(doPrint)fprintf(Stderr,"\n");
298 PolynomialSet h2=g;
299 buchberger(&h2,t);
300 if(doPrint)fprintf(Stderr,"Done computing Gr\"obner basis.\n");
301
302 // AsciiPrinter(Stderr).printPolynomialSet(h2);
303 PolynomialSet wall=initialFormsAssumeMarked(h2,*i);
304 //fprintf(Stderr,"Wall ideal:\n");
305 //AsciiPrinter(Stderr).printPolynomialSet(wall);
306
307 if(containsMonomial(wall))
308 {
309 if(doPrint)fprintf(Stderr,"Initial ideal contains a monomial.\n");
310 Polynomial m(computeTermInIdeal(wall));
311 if(doPrint)fprintf(Stderr,"Done computing term in ideal\n");
312
313 Polynomial temp=m-division(m,h2,LexicographicTermOrder());
314 g.push_back(temp);
315
316 if(doPrint)fprintf(Stderr,"Adding element to basis:\n");
317 if(doPrint)AsciiPrinter(Stderr).printPolynomial(temp);
318 if(doPrint)fprintf(Stderr,"\n");
319
320 *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
321 break;
322 }
323 else
324 {
325 if(doPrint)fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
326 containsNoMonomialCache.insert(*i);
327 }
328 }
329 }
330 if(i==l.end())break;
331 }
332 return g;
333 }
334
335 */
336