1 #include "tropical2.h"
2
3 #include <stdlib.h>
4 #include <iostream>
5
6 #include "buchberger.h"
7 #include "division.h"
8 #include "tropical.h"
9 #include "wallideal.h"
10 #include "dimension.h"
11 #include "halfopencone.h"
12 #include "breadthfirstsearch.h"
13 #include "tropicalbasis.h"
14 #include "tropicalcurve.h"
15 #include "linalg.h"
16
17 #include "field_rationalfunctions2.h"
18
19 #include "timer.h"
20 #include "log.h"
21
22 static Timer tropicalPrincipalIntersectionTimer("Tropical principal intersection",1);
23
24 //////////////////////////////////////////////////////////////////////////
25 //Using stable intersections and rational functions
26
buildH(PolynomialSet const & g,IntegerVector const & isGeneric,bool noHomog=false)27 PolynomialSet buildH(PolynomialSet const &g, IntegerVector const &isGeneric, bool noHomog=false)
28 {
29 int sum=isGeneric.sum();
30 int n=isGeneric.size();
31 PolynomialRing coefRing(g.getRing().getField(),sum);
32 FieldRationalFunctions2 coefField(coefRing);
33 PolynomialRing R(coefField,n-sum);
34
35 PolynomialRing R3=R.withVariablesAppended("H");
36 PolynomialSet ret(R);
37
38 for(PolynomialSet::const_iterator p=g.begin();p!=g.end();p++)
39 {
40 Polynomial q(R);
41 for(TermMap::const_iterator i=p->terms.begin();i!=p->terms.end();i++)
42 {
43 IntegerVector coefv(sum);int IC=0;
44 IntegerVector expv(n-sum);int IE=0;
45 for(int j=0;j<n;j++)
46 if(isGeneric[j])
47 {
48 coefv[IC++]=i->first.exponent[j];
49 }
50 else
51 {
52 expv[IE++]=i->first.exponent[j];
53 }
54 q+=Term(coefField.polynomialToFraction(Polynomial(Term(i->second,Monomial(coefRing,coefv)))),Monomial(R,expv));
55 // Polynomial p(coefRing);
56 // p+=Term(i->m.c,Monomial(coefv));
57 }
58 ret.push_back(q);
59 }
60 if(noHomog)return ret;
61 return ret.homogenization(R3);
62 //We must make sure that the ideal remains homogeneous when turning some variables into parameters.
63 }
64
reconstruct(IntegerVector const & v1,IntegerVector const & isGeneric)65 IntegerVector reconstruct(IntegerVector const &v1, IntegerVector const &isGeneric)
66 {
67 int n=isGeneric.size();
68 IntegerVector v2(n);
69 int I=0;
70 for(int i=0;i<n;i++)if(!isGeneric[i])v2[i]=v1[I++];
71 return v2;
72 }
73
74 /**
75 *
76 */
nonTrivialTropismInner(PolynomialSet const & groebnerBasis)77 IntegerVector nonTrivialTropismInner(PolynomialSet const &groebnerBasis)
78 {
79 PolynomialRing R3=groebnerBasis.getRing().withVariablesAppended("H");
80 PolynomialSet g=saturatedIdeal(groebnerBasis.homogenization(R3));
81 int n=g.getRing().getNumberOfVariables()-1;
82 int h=dimensionOfHomogeneitySpace(g)-1;
83 int d=krullDimension(g)-1;
84
85 // greedy matroid algorithm
86 IntegerVector isGeneric(n);
87 int i=0;
88 while(d>1)
89 {
90 assert(i<n);
91 isGeneric[i]=1;
92
93 bool OK;
94 if(0)
95 { //Using rational functions as coefficients
96 PolynomialSet r=buildH(groebnerBasis,isGeneric);
97 OK=!containsMonomial(r);
98 }
99 else
100 { //Avoiding rational functions as coefficients
101 OK=true;
102 IntegerVectorList M;
103 M.push_back(IntegerVector::allOnes(isGeneric.size())-isGeneric);
104 M.push_back(IntegerVector::allOnes(isGeneric.size()));
105 PolynomialSet G=groebnerBasis;
106 buchberger(&G,MatrixTermOrder(M),true);
107 for(PolynomialSet::const_iterator i=G.begin();i!=G.end();i++)
108 if(!i->isZero())if(i->usedVariables().divides(isGeneric)){
109 OK=false;break;}
110 }
111
112 if(!OK)
113 {
114 isGeneric[i]=0;
115 }
116 else
117 {
118 d--;
119 }
120 i++;
121 }
122 // debug<<"d="<<d<<"\n";
123 if(d==1)
124 {
125 IntegerVectorList A=tropicalCurve(buildH(groebnerBasis,isGeneric,true),true);
126 return reconstruct(A.front(),isGeneric);
127
128 PolynomialSet r=buildH(groebnerBasis,isGeneric);
129 buchberger(&r,LexicographicTermOrder());
130 PolyhedralFan bergmanFan(r.numberOfVariablesInRing());
131 PolynomialSet tropBasis=tropicalBasisOfCurve(r.numberOfVariablesInRing(),r,&bergmanFan,1);
132
133 if(bergmanFan.dimensionOfLinealitySpace()==1)
134 {
135 IntegerVectorList rays=bergmanFan.getRays(2);
136 // At this point the lineality space of the ideal may have increased. A ray is sought outside the original lineality space.
137 assert(!rays.empty());
138 IntegerVector v1=*rays.begin();
139 return reconstruct(v1,isGeneric);
140 }
141 else
142 {
143 assert(bergmanFan.dimensionOfLinealitySpace()==2);
144 IntegerVectorList l=bergmanFan.conesBegin()->generatorsOfLinealitySpace();
145 assert(l.size()==2);
146 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
147 {
148 if(!groebnerBasis.isHomogeneous(reconstruct(*i,isGeneric)))
149 return reconstruct(*i,isGeneric);
150 }
151 pout<<l;
152 pout<<groebnerBasis;
153 assert(0);
154 }
155 }
156 return IntegerVector(n);
157 }
158
nonTrivialTropism(PolynomialSet const & groebnerBasis)159 IntegerVector nonTrivialTropism(PolynomialSet const &groebnerBasis)
160 {
161 list<int> toKeep=groebnerBasis.multiDeHomogenizationToKeep();
162 PolynomialSet ideal2=groebnerBasis.multiDeHomogenization();
163 IntegerVector v=nonTrivialTropismInner(ideal2);
164 IntegerVector ret(groebnerBasis.numberOfVariablesInRing());
165 int I=0;
166 for(list<int>::const_iterator i=toKeep.begin();i!=toKeep.end();i++,I++)
167 ret[*i]=v[I];
168 return ret;
169 }
170
startingConeError()171 static void startingConeError()
172 {
173 fprintf(Stderr,"UNABLE TO COMPUTE STARTING CONE.\n");
174 fprintf(Stderr,"THE STARTING CONE ALGORITHM IN GFAN IS BASED ON HEURISTICS WHICH HAVE FAILED ON THIS EXAMPLE.\n");
175 assert(0);
176 }
177
initialIdeal(PolynomialSet const & g,IntegerVector const & weight)178 PolynomialSet initialIdeal(PolynomialSet const &g, IntegerVector const &weight)
179 //Assume homogeneous
180 {
181 PolynomialSet ret=g;
182 WeightReverseLexicographicTermOrder T(weight);
183 buchberger(&ret,T);
184 return initialForms(ret,weight);
185 }
186
initialIdealNonHomogeneous(PolynomialSet const & g,IntegerVector const & weight,bool allowSaturation)187 PolynomialSet initialIdealNonHomogeneous(PolynomialSet const &g, IntegerVector const &weight, bool allowSaturation)
188 {
189 // We use Theorem 4 page 379 of [Cox, Little, O'Shea] (book 1) to compute the homogenization
190 IntegerVector w(IntegerVector::allOnes(g.getRing().getNumberOfVariables()));
191 WeightReverseLexicographicTermOrder T(w);
192 PolynomialSet g2=g;
193 buchberger(&g2,T,allowSaturation);
194 PolynomialRing R2(g.getRing().getField(),g.getRing().getNumberOfVariables()+1);
195 PolynomialSet g3=g2.homogenization(R2,&w);
196 // We now compute the initial ideal using Proposition 5.2.3 in [Jensen] (PhD thesis).
197 IntegerVector weight2=concatenation(weight,IntegerVector(1));
198 WeightTermOrder/*WeightReverseLexicographicTermOrder*/ T2(weight2);// The proposition assumes that weight2 is tie-broken with a term order, but that is probably not necessary. However, the specifications of _this_ routine says that a GB must be returned.
199 buchberger(&g3,T2,allowSaturation);
200 PolynomialSet g4=initialFormsAssumeMarked(g3,weight2);
201 return g4.deHomogenization();
202 }
203
204
initialFormAssumeMarked(Polynomial const & p,IntegerVector const & weight)205 Polynomial initialFormAssumeMarked(Polynomial const &p, IntegerVector const &weight)
206 {
207 Polynomial r(p.getRing());
208 IntegerVector markedExponent=p.getMarked().m.exponent;
209
210 for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
211 {
212 IntegerVector dif=markedExponent-i->first.exponent;
213
214 if(dot(dif,weight)==0)
215 r+=Polynomial(Term(i->second,i->first));
216 }
217 r.mark(Monomial(p.getRing(),markedExponent));
218
219 return r;
220 }
221
222
initialFormsAssumeMarked(PolynomialSet const & groebnerBasis,IntegerVector const & weight)223 PolynomialSet initialFormsAssumeMarked(PolynomialSet const &groebnerBasis, IntegerVector const &weight)
224 {
225 PolynomialRing theRing=groebnerBasis.getRing();
226 PolynomialSet r(theRing);
227
228 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
229 {
230 r.push_back(initialFormAssumeMarked(*i,weight));
231 }
232 return r;
233 }
234
235
initialForm(Polynomial const & p,IntegerVector const & weight)236 Polynomial initialForm(Polynomial const &p, IntegerVector const &weight)
237 {
238 if(p.isZero())return p;
239 int64 a=dotLong(p.terms.begin()->first.exponent,weight);
240 for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
241 {
242 int64 b=dotLong(i->first.exponent,weight);
243 if(b>a)a=b;
244 }
245
246 Polynomial r(p.getRing());
247 bool ismarked=p.isMarked();
248 IntegerVector markedExponent;
249 if(ismarked)markedExponent=p.getMarked().m.exponent;
250
251 bool markedFound=false;
252
253 for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
254 {
255 if(dotLong(i->first.exponent,weight)==a)
256 {
257 r+=Polynomial(Term(i->second,i->first));
258 if(ismarked)if((markedExponent-(i->first.exponent)).isZero())markedFound=true;
259 }
260 }
261 if(markedFound)
262 r.mark(Monomial(p.getRing(),markedExponent));
263
264 return r;
265 }
266
267
initialForms(PolynomialSet const & groebnerBasis,IntegerVector const & weight)268 PolynomialSet initialForms(PolynomialSet const &groebnerBasis, IntegerVector const &weight)
269 {
270 PolynomialRing theRing=groebnerBasis.getRing();
271 PolynomialSet r(theRing);
272 if(theRing.getNumberOfVariables()!=weight.size())
273 {
274 cerr << "Error: Number of varaibles in polynomial ring "<<theRing.getNumberOfVariables()<< " length of weight vector " << weight.size() <<endl;
275 assert(0);
276 }
277
278 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
279 {
280 r.push_back(initialForm(*i,weight));
281 }
282 return r;
283 }
284
285
286
tropicalPrincipalIntersection(int n,PolynomialSet const & g,int linealitySpaceDimension)287 PolyhedralFan tropicalPrincipalIntersection(int n, PolynomialSet const &g, int linealitySpaceDimension)
288 {
289 //return tropicalHyperSurfaceIntersection(n, g);////////////////////////////////////////
290 gfan_log2 fprintf(Stderr,"Intersecting\n");
291 log3 AsciiPrinter(Stderr).printPolynomialSet(g);
292
293 TimerScope ts(&tropicalPrincipalIntersectionTimer);
294 PolyhedralFan ret=PolyhedralFan::fullSpace(n);
295
296 for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
297 {
298 ret=refinement(ret,PolyhedralFan::bergmanOfPrincipalIdeal(*i),linealitySpaceDimension,true);
299 }
300 gfan_log2 fprintf(Stderr,"Done intersecting\n");
301 return ret;
302 }
303
304
305
checkList(IntegerVectorList const & l,PolynomialSet const & groebnerBasis,PolynomialSet * fullNeighbourBasis,int h,bool & result,bool onlyCheckRays)306 static PolynomialSet checkList(IntegerVectorList const &l, PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, int h, bool &result, bool onlyCheckRays)
307 {
308 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
309 {
310 WeightReverseLexicographicTermOrder t(*i);
311 gfan_log2 fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
312 gfan_log2 AsciiPrinter(Stderr).printVector(*i);
313 gfan_log2 fprintf(Stderr,"\n");
314 PolynomialSet h2=groebnerBasis;
315 buchberger(&h2,t);
316 gfan_log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
317
318 log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
319 PolynomialSet wall=initialFormsAssumeMarked(h2,*i);
320
321 log3 AsciiPrinter(Stderr).printString("Initial ideal:\n");
322 log3 AsciiPrinter(Stderr).printPolynomialSet(wall);
323
324 int hdim2=dimensionOfHomogeneitySpace(wall);
325 if(hdim2>h)
326 {
327 if(!containsMonomial(wall))
328 {
329 log1 fprintf(Stderr,"Iterating recursively.\n");
330 //PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,0);
331 PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,fullNeighbourBasis,onlyCheckRays);
332
333 if(fullNeighbourBasis)
334 {
335 //*fullNeighbourBasis=liftBasis(initialIdeal,h2);
336 *fullNeighbourBasis=liftBasis(*fullNeighbourBasis,h2);
337 }
338
339
340 result=true;
341 return initialIdeal;
342 }
343 }
344 }
345 result=false;
346 return groebnerBasis;
347 }
348
349
350
perturbedRelativeInteriorTropism(PolynomialSet const & groebnerBasis)351 IntegerVectorList perturbedRelativeInteriorTropism(PolynomialSet const &groebnerBasis)
352 {
353 int n=groebnerBasis.getRing().getNumberOfVariables();
354 // If T(I) is just the origin, we return the empty list
355 PolynomialRing R(groebnerBasis.getRing().getField(),n+1);
356 PolynomialSet g=groebnerBasis.homogenization(R);
357 saturatedIdeal(g);
358 if(krullDimension(g)==1)return IntegerVectorList();
359
360 // If the homogeneity space is non-trivial find a basis of it and reduce dimension
361 PolyhedralCone hom=homogeneitySpace(g);
362 if(hom.dimension()>1)
363 {
364 FieldMatrix m=integerMatrixToFieldMatrix(rowsToIntegerMatrix(hom.getEquations(),n+1),Q);
365 m.reduce();
366 list<int> toKeep=m.pivotColumns();
367 PolynomialRing R2=PolynomialRing(g.getRing().getField(),toKeep.size());
368
369 PolynomialSet i2=g.embeddedInto(R2,&toKeep);
370
371 IntegerVectorList rl=perturbedRelativeInteriorTropism(i2);
372
373 IntegerVectorList eq;eq.push_back(IntegerVector::standardVector(n+1,n));
374 PolyhedralCone C(IntegerVectorList(),eq,n+1);C=intersection(C,hom);C.canonicalize();
375 IntegerVectorList ret=C.projection(n).generatorsOfLinealitySpace();
376 for(IntegerVectorList::const_iterator j=rl.begin();j!=rl.end();j++)
377 {
378 IntegerVector temp(n+1);
379 int I=0;
380 for(list<int>::const_iterator i=toKeep.begin();i!=toKeep.end();i++,I++)
381 temp[*i]=(*j)[I];
382 ret.push_back(temp.subvector(0,n)-temp[n]*IntegerVector::allOnes(n));
383 }
384 return ret;
385 }
386
387 IntegerVector v=nonTrivialTropism(groebnerBasis);
388 IntegerVectorList ret=perturbedRelativeInteriorTropism(initialIdealNonHomogeneous(groebnerBasis,v,true)); //At this point it would make sense to "dehomogenize w.r.t. v", so that v is not repeated in the recursion. We could also simply choose not to push v.
389 ret.push_front(v);
390
391 return ret;
392
393 /* saturate
394 if dim=0 return empty;
395 compute homogeneity space
396 if homog>0
397 reduce dimension of ring/ideal
398 call recursively
399 lift list
400 return
401 else
402 return nontrivialtropism
403 */
404 }
405
406
guessInitialIdealWithoutMonomial(PolynomialSet const & groebnerBasis,PolynomialSet * fullNeighbourBasis,bool onlyCheckRays)407 PolynomialSet guessInitialIdealWithoutMonomial(PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, bool onlyCheckRays) //ideal must be homogeneous
408 // fullNeighbourBasis is set to a Groebner basis of the full ideal. The returned basis and fullNeighbourBasis have at least one termorder in common
409 {
410 if(1)
411 {
412 IntegerVectorList w=perturbedRelativeInteriorTropism(groebnerBasis);
413 MatrixTermOrder T(w);
414 PolynomialSet g=groebnerBasis;
415 buchberger(&g,T);//assuming that g is homogeneous
416 PolynomialSet ig=g;
417 for(IntegerVectorList::const_iterator i=w.begin();i!=w.end();i++)ig=initialForms(ig,*i);
418 assert(fullNeighbourBasis);
419 *fullNeighbourBasis=g;
420 return ig;
421 }
422 // log0 fprintf(Stderr,"A\n");
423 assert(groebnerBasis.isValid());
424 // log0 fprintf(Stderr,"B\n");
425 if(fullNeighbourBasis)
426 {
427 assert(fullNeighbourBasis->isValid());
428 }
429 // log0 fprintf(Stderr,"C\n");
430
431 int n=groebnerBasis.numberOfVariablesInRing();
432 // log0 fprintf(Stderr,"D\n");
433 int h=dimensionOfHomogeneitySpace(groebnerBasis);
434 // log0 fprintf(Stderr,"E\n");
435 int d=krullDimension(groebnerBasis);
436 // log0 fprintf(Stderr,"F\n");
437
438 if(d==h)
439 {
440 if(fullNeighbourBasis)*fullNeighbourBasis=groebnerBasis;
441 return groebnerBasis;
442 }
443
444 #if 0
445 // stable intersections/rational functions
446 {
447 IntegerVectorList toCheck;
448 // pout<<groebnerBasis;
449
450
451 toCheck.push_back(nonTrivialTropism(groebnerBasis));
452 bool result;
453 pout<<toCheck<<n<<h<<d<<"\n";
454 PolynomialSet r=checkList(toCheck,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
455 if(result)return r;
456 }
457 assert(0);
458 // TODO: When starting cone using stable intersections works, remove all dead code.
459 #endif
460
461 {
462 //gfan_log2
463 fprintf(Stderr,"Computing extreme rays.\n");
464 //IntegerVectorList a;
465 PolyhedralCone p=coneFromMarkedBasis(groebnerBasis);
466 //PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
467 IntegerVectorList extreme=p.extremeRays();
468 gfan_log2 fprintf(Stderr,"Extreme rays of Groebner cone:\n");
469 gfan_log2 AsciiPrinter(Stderr).printVectorList(extreme);
470
471 bool result;
472 PolynomialSet r=checkList(extreme,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
473 if(result)return r;
474 }
475 if(onlyCheckRays)startingConeError();
476 PolyhedralFan f=PolyhedralFan::fullSpace(n);
477 /* for(int i=0;i<d-1;i++)
478 {
479 IntegerVector v(n);
480 for(int j=0;j<n;j++)v[j]=rand()&1;
481 IntegerVectorList a,b;
482 b.push_back(v);
483 PolyhedralFan F(n);
484 F.insert(PolyhedralCone(a,b));
485 f=refinement(f,F);
486 }
487 AsciiPrinter P(Stderr);
488 f.print(&P);
489 */
490 int hypersurfacesToGo=groebnerBasis.size();
491 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
492 {
493 fprintf(Stderr,"Hypersurfaces to go:%i\n",hypersurfacesToGo--);
494 fprintf(Stderr,"Max dimension: %i\n",f.getMaxDimension());
495 debug<<"*i="<<*i<<"\n";
496 f=refinement(f,PolyhedralFan::bergmanOfPrincipalIdeal(*i));
497 f.removeAllExcept(3);
498
499 IntegerVectorList l=f.getRelativeInteriorPoints();
500
501 bool result;
502 PolynomialSet r=checkList(l,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
503 if(result)return r;
504 }
505 startingConeError();
506 return groebnerBasis;
507 }
508
checkListStably(IntegerVectorList const & l,PolynomialSet const & groebnerBasis,PolynomialSet * fullNeighbourBasis,int h,bool & result,bool onlyCheckRays)509 static PolynomialSet checkListStably(IntegerVectorList const &l, PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, int h, bool &result, bool onlyCheckRays)
510 {
511 debug<< "Checklist called on"<<groebnerBasis;
512 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
513 {
514 WeightReverseLexicographicTermOrder t(*i);
515 gfan_log2 fprintf(Stderr,"Taking initial forms with respect to:");
516 gfan_log2 AsciiPrinter(Stderr).printVector(*i);
517 gfan_log2 fprintf(Stderr,"\n");
518 PolynomialSet h2=groebnerBasis;
519 gfan_log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
520
521 log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
522 PolynomialSet wall=initialForms(h2,*i);
523
524 log3 AsciiPrinter(Stderr).printString("Initial ideal:\n");
525 log3 AsciiPrinter(Stderr).printPolynomialSet(wall);
526
527 int hdim2=dimensionOfHomogeneitySpace(wall);
528 if(hdim2>h)
529 {
530 if(nonEmptyStableIntersection(wall))
531 {
532 log1 fprintf(Stderr,"Iterating recursively.\n");
533 //PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,0);
534 PolynomialSet initialIdeal=guessInitialIdealWithoutMonomialStably(wall,fullNeighbourBasis,onlyCheckRays);
535
536 if(fullNeighbourBasis)
537 {
538 //*fullNeighbourBasis=liftBasis(initialIdeal,h2);
539 // *fullNeighbourBasis=liftBasis(*fullNeighbourBasis,h2);
540 *fullNeighbourBasis=groebnerBasis;
541 fullNeighbourBasis->copyMarkings(initialIdeal);
542 }
543
544
545 result=true;
546 return initialIdeal;
547 }
548 }
549 }
550 result=false;
551 return groebnerBasis;
552 }
553
guessInitialIdealWithoutMonomialStably(PolynomialSet const & groebnerBasis,PolynomialSet * fullNeighbourBasis,bool onlyCheckRays)554 PolynomialSet guessInitialIdealWithoutMonomialStably(PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, bool onlyCheckRays) //ideal must be homogeneous
555 // fullNeighbourBasis is set to a Groebner basis of the full ideal. The returned basis and fullNeighbourBasis have at least one termorder in common
556 {
557 int n=groebnerBasis.numberOfVariablesInRing();
558 int h=dimensionOfHomogeneitySpace(groebnerBasis);
559 int d=n-groebnerBasis.size();//krullDimension(groebnerBasis);
560
561 debug<</*"d"<<d<<*/"h"<<h<<"n"<<n<<"\n";
562
563
564 if(d==h)
565 {
566 if(fullNeighbourBasis)*fullNeighbourBasis=groebnerBasis;
567 return groebnerBasis;
568 }
569
570 {
571 gfan_log2 fprintf(Stderr,"Computing extreme rays.\n");
572 //IntegerVectorList a;
573 PolyhedralCone p=coneFromMarkedBasis(groebnerBasis);
574 //PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
575 IntegerVectorList extreme=p.extremeRays();
576 gfan_log2 fprintf(Stderr,"Extreme rays of Groebner cone:\n");
577 gfan_log2 AsciiPrinter(Stderr).printVectorList(extreme);
578
579 bool result;
580 PolynomialSet r=checkListStably(extreme,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
581 if(result)return r;
582 }
583 if(onlyCheckRays)startingConeError();
584
585 PolyhedralFan f=PolyhedralFan::fullSpace(n);
586
587 int hypersurfacesToGo=groebnerBasis.size();
588 for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
589 {
590 fprintf(Stderr,"Hypersurfaces to go:%i\n",hypersurfacesToGo--);
591 fprintf(Stderr,"Max dimension: %i\n",f.getMaxDimension());
592 f=refinement(f,PolyhedralFan::bergmanOfPrincipalIdeal(*i));
593 f.removeAllExcept(3);
594
595 IntegerVectorList l=f.getRelativeInteriorPoints();
596
597 bool result;
598 PolynomialSet r=checkListStably(l,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
599 if(result)return r;
600 }
601 startingConeError();
602 return groebnerBasis;
603 }
604