14
18
typedef set<IntegerVector> IntegerVectorSet;
20
static Polynomial cleverSaturation(Polynomial const &p, IntegerVector const &w)
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();
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());
37
Polynomial all=theRing.monomialFromExponentVector(s*IntegerVector::allOnes(theRing.getNumberOfVariables())-gcd);
38
Polynomial g= division(all*p,reducer,LexicographicTermOrder());
40
debug<<"NEW"<<g<<"\n";
46
static void initialSaturatingBuchberger(PolynomialSet *g, TermOrder const &termOrder_, IntegerVector const &w_)
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();
53
MatrixTermOrder termOrder(temp);
54
PolynomialRing theRing=g->getRing().withVariablesAppended("Z");
55
PolynomialSet sPolynomials(theRing);
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
60
sPolynomials.push_back(theRing.monomialFromExponentVector(IntegerVector::allOnes(theRing.getNumberOfVariables()))-theRing.one());
62
sPolynomials.saturate();
63
sPolynomials.markAndScale(termOrder);
65
*g=PolynomialSet(theRing);
67
while(!sPolynomials.empty())
69
Polynomial p=*sPolynomials.begin();
70
sPolynomials.pop_front();
72
p=division(p,*g,termOrder);
75
p=cleverSaturation(p,w);
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!
82
if(!relativelyPrime(i->getMarked().m.exponent,p.getMarked().m.exponent))
84
Polynomial s=sPolynomial(*i,p);
85
s.mark(termOrder); // with respect to some termorder
86
s.scaleMarkedCoefficientToOne();
87
sPolynomials.push_back(s);
94
// if((t&31)==0)fprintf(Stderr," gsize %i spolys:%i\n",g->size(),sPolynomials.size());
99
autoReduce(g,termOrder);
16
104
PolynomialSet tropicalBasisOfCurve(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension) //Assuming g is homogeneous
18
int homog=linealitySpaceDimension;
108
* TODO: Introduce the notion of a tropical prebasis:
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.
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.
118
// bool prebasis=true;
119
// debug<<"PREBASIS="<<prebasis<<"\n";
120
log2 debug<<"TropicalBasis begin\n";
122
int homog=linealitySpaceDimension;
19
123
assert(homog>0 || n==0);
20
124
TimerScope ts(&iterativeTropicalBasisTimer);
21
125
PolyhedralFan f(n);
22
126
if(!intersectionFan)intersectionFan=&f;
24
*intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);
128
// *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);
129
// log1 fprintf(Stderr,"WARINING USING EXPERIMENTAL TROPICAL HYPERSURFACE INTERSECTION ROUTINE!!\n");
130
*intersectionFan=tropicalHyperSurfaceIntersectionClosed(n, g);
26
132
IntegerVectorSet containsNoMonomialCache;
30
136
PolyhedralFan::coneIterator i;
138
// {AsciiPrinter P(Stderr);intersectionFan->printWithIndices(&P);}
140
// {AsciiPrinter P(Stderr);intersectionFan->printWithIndices(&P);}
31
141
for(i=intersectionFan->conesBegin();i!=intersectionFan->conesEnd();i++)
143
// log1 cerr<<"!@#$";
33
144
IntegerVector relativeInteriorPoint=i->getRelativeInteriorPoint();
145
// log1 cerr<<"1234/n";
34
147
if(containsNoMonomialCache.count(relativeInteriorPoint)>0)
36
149
log2 fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
155
if(containsMonomial(initialForms(g,relativeInteriorPoint)))
157
intersectionFan->insertFacetsOfCone(*i);
158
intersectionFan->remove(*i);
159
debug<<"LOWERING DIMENSION OF CONE\n";//TODO: checking cones in order of dimension could avoid this.
40
163
WeightReverseLexicographicTermOrder t(relativeInteriorPoint);
41
164
log2 fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
42
165
log2 AsciiPrinter(Stderr).printVector(relativeInteriorPoint);
43
166
log2 fprintf(Stderr,"\n");
44
167
PolynomialSet h2=g;
170
/* {Adjust these lines somehow to enable the saturating buchberger.
172
initialSaturatingBuchberger(&h2, t, relativeInteriorPoint);
173
//buchberger(&h2,t,true);
174
debug<<"SATURATING BUCHBERGER DONE\n";
176
// buchberger(&h2,t);
177
h2=GE_groebnerBasis(h2,t,true);
46
179
log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");
48
log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
182
// log3 AsciiPrinter(Stderr).printPolynomialSet(h2);
50
184
PolynomialSet wall=initialFormsAssumeMarked(h2,relativeInteriorPoint);
52
186
if(containsMonomial(wall))
54
188
log2 fprintf(Stderr,"Initial ideal contains a monomial.\n");
55
189
Polynomial m(computeTermInIdeal(wall));
56
190
log2 fprintf(Stderr,"Done computing term in ideal\n");
58
192
Polynomial temp=m-division(m,h2,LexicographicTermOrder());
59
193
g.push_back(temp);
61
195
log2 fprintf(Stderr,"Adding element to basis:\n");
62
196
log2 AsciiPrinter(Stderr).printPolynomial(temp);
63
197
log2 fprintf(Stderr,"\n");
65
199
*intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
added
removed
tropicalbasis.cpp
