← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to integergb.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
integergb.cpp
 
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
 
 
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
 
 
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
 
 
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
 
 
128
Polynomial spol(Polynomial const &g1, Polynomial const &g2)
 
129
{
 
130
  return sgpol(g1,g2,true);
 
131
}
 
132
 
 
133
 
 
134
Polynomial gpol(Polynomial const &g1, Polynomial const &g2)
 
135
{
 
136
  return sgpol(g1,g2,false);
 
137
}
 
138
 
 
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
 
 
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;
 
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
 
 
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
 
 
273
IntegerGroebnerFanTraverser::IntegerGroebnerFanTraverser(PolynomialSet const &generators):
 
274
  n(generators.getRing().getNumberOfVariables()),
 
275
  groebnerBasis(generators)
 
276
{
 
277
//  LexicographicTermOrder tieBreaker;
 
278
//  zBuchberger(groebnerBasis,tieBreaker);
 
279
//debug<<"WE ASSUME THAT WE ALREADY HAVE A REDUCED GB\n";
 
280
 
 
281
  updatePolyhedralCone();
 
282
}
 
283
 
 
284
void IntegerGroebnerFanTraverser::changeCone(IntegerVector const &ridgeVector, IntegerVector const &rayVector)
 
285
{
 
286
  IntegerVectorList tOld;
 
287
  tOld.push_back(ridgeVector);
 
288
  tOld.push_back(-rayVector);
 
289
  MatrixTermOrder TOld(tOld);
 
290
 
 
291
  IntegerVectorList t;
 
292
  t.push_back(ridgeVector);
 
293
  t.push_back(rayVector);
 
294
  MatrixTermOrder T(t);
 
295
 
 
296
//debug<<"FLIP!\n";
 
297
//debug<<"OLDGB\n"<<groebnerBasis;
 
298
PolynomialSet g=initialFormsAssumeMarked(groebnerBasis,ridgeVector);
 
299
//debug<<"OLD INITIAL GB\n"<<g;
 
300
zBuchberger(g,T);
 
301
//debug<<"NEW INITIAL GB\n"<<g;
 
302
 
 
303
 
 
304
 
 
305
PolynomialSet liftedBasis(groebnerBasis.getRing());
 
306
for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
 
307
  {
 
308
    Polynomial lt=*i;//Polynomial(i->getMarked());
 
309
    Polynomial r=eDivision(lt, groebnerBasis,TOld);
 
310
    Polynomial rIn=eDivision(lt,initialFormsAssumeMarked(groebnerBasis,ridgeVector),TOld);
 
311
    liftedBasis.push_back(lt-r);
 
312
 
 
313
//   debug<<"--------------------\n";
 
314
//    debug<<"lt:"<<lt<<"\n";
 
315
//    debug<<"r:"<<r<<"\n";
 
316
//    debug<<"rIn:"<<rIn<<"\n";
 
317
//    debug<<"--------------------\n";
 
318
  }
 
319
liftedBasis.mark(T);
 
320
//debug<<"to be autoreduced:"<<liftedBasis;
 
321
//debug<<"autoreducing......\n";
 
322
liftedBasis.mark(T);
 
323
zAutoReduce(&liftedBasis,T);
 
324
PolynomialSet newInitialForms=initialForms(liftedBasis,ridgeVector);
 
325
//debug<<"LIFTED"<<liftedBasis;
 
326
//debug<<"NEWINITIALFORMS"<<newInitialForms;
 
327
 
 
328
/*zBuchberger(groebnerBasis,T);
 
329
 
 
330
{
 
331
  groebnerBasis.sort_();
 
332
  liftedBasis.sort_();
 
333
  bool areEqual=(groebnerBasis.size()==liftedBasis.size());
 
334
  if(areEqual)
 
335
    {
 
336
      PolynomialSet::const_iterator i=groebnerBasis.begin();
 
337
      for(PolynomialSet::const_iterator j=liftedBasis.begin();j!=liftedBasis.end();j++,i++)
 
338
        areEqual&=(*i-*j).isZero();
 
339
    }
 
340
  if(!areEqual)
 
341
    {
 
342
      debug<<groebnerBasis<<liftedBasis;
 
343
      assert(0);
 
344
    }
 
345
}*/
 
346
 
 
347
groebnerBasis=liftedBasis;
 
348
groebnerBasis.mark(T);
 
349
 
 
350
 
 
351
//debug<<"NEW BASIS"<<groebnerBasis;
 
352
 
 
353
  updatePolyhedralCone();
 
354
}
 
355
 
 
356
IntegerVectorList IntegerGroebnerFanTraverser::link(IntegerVector const &ridgeVector)
 
357
{
 
358
  IntegerVectorList ret;
 
359
  IntegerVector v=theCone.link(ridgeVector).getUniquePoint();
 
360
  ret.push_back(v);
 
361
 
 
362
  PolyhedralCone temp=intersection(PolyhedralCone::positiveOrthant(n),theCone.faceContaining(ridgeVector));
 
363
  IntegerVector temp2=temp.getRelativeInteriorPoint();
 
364
  if(temp2.min()>0)
 
365
    {
 
366
      ret.push_back(-v);
 
367
    }
 
368
//  debug<<"LINK"<<ret;
 
369
  return ret;
 
370
}
 
371
 
 
372
 
 
373
PolyhedralCone & IntegerGroebnerFanTraverser::refToPolyhedralCone()
 
374
{
 
375
  return theCone;
 
376
}