gitfan.cc
Go to the documentation of this file.
1 /***************************************************************
2  *
3  * File: gitfan.cc
4  * Purpose: Computationally intensive procedures for gitfan.lib,
5  * outsourced to improve the performance.
6  * Authors: Janko Boehm boehm@mathematik.uni-kl.de
7  * Simon Keicher keicher@mail.mathematik.uni-tuebingen.de
8  * Yue Ren ren@mathematik.uni-kl.de
9  *
10  ***************************************************************/
11 
12 #include <kernel/mod2.h>
13 
14 #if HAVE_GFANLIB
15 
16 #include <callgfanlib_conversion.h>
17 #include <bbcone.h>
18 #include <bbfan.h>
19 #include <gitfan.h>
20 
21 #include <Singular/ipid.h>
22 #include <Singular/lists.h>
23 #include <Singular/ipshell.h>
24 
25 #include <coeffs/bigintmat.h>
26 
27 
28 namespace gitfan
29 {
30 
32  eta(gfan::ZCone()),
33  interiorPoint(gfan::ZVector()),
34  facetNormal(gfan::ZVector())
35  {
36  }
37 
39  eta(f.eta),
40  interiorPoint(f.interiorPoint),
41  facetNormal(f.facetNormal)
42  {
43 #ifndef SING_NDEBUG
44  gfan::ZCone c = f.eta;
45  gfan::ZVector v = f.interiorPoint;
46  gfan::ZVector w = f.facetNormal;
47  assume(c.ambientDimension() == (int)v.size());
48  assume(c.ambientDimension() == (int)w.size());
49  assume(c.contains(v));
50  assume(!c.contains(w));
51 #endif
52  }
53 
54  facet::facet(const gfan::ZCone &c, const gfan::ZVector &v, const gfan::ZVector &w):
55  eta(c),
56  interiorPoint(v),
57  facetNormal(w)
58  {
59 #ifndef SING_NDEBUG
60  assume(c.ambientDimension() == (int)v.size());
61  assume(c.ambientDimension() == (int)w.size());
62  assume(c.contains(v));
63  assume(!c.contains(w));
64 #endif
65  }
66 
68  {
69 #ifndef SING_NDEBUG
70  gfan::ZCone c = this->eta;
71  gfan::ZVector v = this->interiorPoint;
72  gfan::ZVector w = this->facetNormal;
73  assume(c.ambientDimension() == (int)v.size());
74  assume(c.ambientDimension() == (int)w.size());
75  assume(c.contains(v));
76  assume(!c.contains(w));
77 #endif
78  }
79 
80  void mergeFacets(facets &F, const facets &newFacets)
81  {
82  std::pair<facets::iterator,bool> check(newFacets.begin(),false);
83  for(facets::iterator p=newFacets.begin(); p!=newFacets.end(); p++)
84  {
85  check = F.insert(*p);
86  if(!check.second)
87  F.erase(check.first);
88  }
89  }
90 
91 }
92 
93 
94 static gfan::ZCone subcone(const lists &cones, const gfan::ZVector &point)
95 {
96  gfan::ZCone sigma = gfan::ZCone(gfan::ZMatrix(1,point.size()), gfan::ZMatrix(1,point.size()));
97  gfan::ZCone* zc;
98  for (int i=0; i<=cones->nr; i++)
99  {
100  zc = (gfan::ZCone*) cones->m[i].Data();
101  if (zc->contains(point))
102  sigma = gfan::intersection(sigma,*zc);
103  }
104  return(sigma);
105 }
106 
107 static gitfan::facets interiorFacets(const gfan::ZCone &zc, const gfan::ZCone &bound)
108 {
109  gfan::ZMatrix inequalities = zc.getFacets();
110  gfan::ZMatrix equations = zc.getImpliedEquations();
111  int r = inequalities.getHeight();
112  int c = inequalities.getWidth();
113  gitfan::facets F;
114  if (r*c == 0)
115  /***
116  * this is the trivial case where either we are in a zerodimensional ambient space,
117  * or the cone has no facets.
118  **/
119  return F;
120 
121  // int index = 0;
122  /* next we iterate over each of the r facets, build the respective cone and add it to the list */
123  /* this is the i=0 case */
124  gfan::ZMatrix newInequalities = inequalities.submatrix(1,0,r,c);
125  gfan::ZMatrix newEquations = equations;
126  newEquations.appendRow(inequalities[0]);
127  gfan::ZCone eta = gfan::ZCone(newInequalities,newEquations);
128  eta.canonicalize();
129  gfan::ZVector v = eta.getRelativeInteriorPoint();
130  gfan::ZVector w = inequalities[0];
131 
132  if (bound.containsRelatively(v))
133  F.insert(gitfan::facet(eta,v,w));
134 
135  /* these are the cases i=1,...,r-2 */
136  for (int i=1; i<r-1; i++)
137  {
138  newInequalities = inequalities.submatrix(0,0,i,c);
139  newInequalities.append(inequalities.submatrix(i+1,0,r,c));
140  newEquations = equations;
141  newEquations.appendRow(inequalities[i]);
142  eta = gfan::ZCone(newInequalities,newEquations);
143  eta.canonicalize();
144  v = eta.getRelativeInteriorPoint();
145  w = inequalities[i];
146  if (bound.containsRelatively(v))
147  F.insert(gitfan::facet(eta,v,w));
148  }
149 
150  /* this is the i=r-1 case */
151  newInequalities = inequalities.submatrix(0,0,r-1,c);
152  newEquations = equations;
153  newEquations.appendRow(inequalities[r-1]);
154  eta = gfan::ZCone(newInequalities,newEquations);
155  eta.canonicalize();
156 
157  v = eta.getRelativeInteriorPoint();
158  w = inequalities[r-1];
159  if (bound.containsRelatively(v))
160  F.insert(gitfan::facet(eta,v,w));
161 
162  return F;
163 }
164 
166 {
167  leftv u=args;
168  if ((u != NULL) && (u->Typ() == LIST_CMD))
169  {
170  leftv v=u->next;
171  if ((v != NULL) && (v->Typ() == BIGINTMAT_CMD))
172  {
173  lists cones = (lists) u->Data();
174  bigintmat* bim = (bigintmat*) v->Data();
175  gfan::ZMatrix* zm = bigintmatToZMatrix(bim->transpose());
176  gfan::ZCone support = gfan::ZCone::givenByRays(*zm, gfan::ZMatrix(0, zm->getWidth()));
177  delete zm;
178 
179  /***
180  * Randomly compute a first full-dimensional cone and insert it into the fan.
181  * Compute a list of facets and relative interior points.
182  * The relative interior points are unique, assuming the cone is stored in canonical form,
183  * which is the case in our algorithm, as we supply no redundant inequalities.
184  * Hence we can decide whether a facet need to be traversed by crosschecking
185  * its relative interior point with this list.
186  **/
187  gfan::ZCone lambda; gfan::ZVector point;
188  do
189  {
190  point = randomPoint(&support);
191  lambda = subcone(cones, point);
192  }
193  while (lambda.dimension() < lambda.ambientDimension());
194  int iterationNumber = 1;
195  std::cout << "cones found: " << iterationNumber++ << std::endl;
196 
197  lambda.canonicalize();
198  gfan::ZFan* Sigma = new gfan::ZFan(lambda.ambientDimension());
199  Sigma->insert(lambda);
200  gitfan::facets F = interiorFacets(lambda, support);
201  if (F.empty())
202  {
203  res->rtyp = fanID;
204  res->data = (void*) Sigma;
205  return FALSE;
206  }
207  int mu = 1024;
208 
210  gfan::ZCone eta;
211  gfan::ZVector interiorPoint;
212  gfan::ZVector facetNormal;
213  gitfan::facets newFacets;
214  while (!F.empty())
215  {
216  /***
217  * Extract a facet to traverse and its relative interior point.
218  **/
219  f = *(F.begin());
220  eta = f.getEta();
221  interiorPoint = f.getInteriorPoint();
222  facetNormal = f.getFacetNormal();
223 
224  /***
225  * construct a point, which lies on the other side of the facet.
226  * make sure it lies in the known support of our fan
227  * and that the cone around the point is maximal, containing eta.
228  **/
229  point = mu * interiorPoint - facetNormal;
230  while (!support.containsRelatively(point))
231  {
232  mu = mu * 16;
233  point = mu * interiorPoint - facetNormal;
234  }
235 
236  lambda = subcone(cones,point);
237  while ((lambda.dimension() < lambda.ambientDimension()) && !(lambda.contains(interiorPoint)))
238  {
239  mu = mu * 16;
240  point = mu * interiorPoint - facetNormal;
241  lambda = subcone(cones,point);
242  }
243  std::cout << "cones found: " << iterationNumber++ << std::endl;
244 
245  /***
246  * insert lambda into Sigma, and create a list of facets of lambda.
247  * merge the two lists of facets
248  **/
249  lambda.canonicalize();
250  Sigma->insert(lambda);
251  newFacets = interiorFacets(lambda, support);
252  mergeFacets(F,newFacets);
253  newFacets.clear();
254  }
255  res->rtyp = fanID;
256  res->data = (void*) Sigma;
257  return FALSE;
258  }
259  }
260  WerrorS("refineCones: unexpected parameters");
261  return TRUE;
262 }
263 
264 
265 static int binomial(int n, int k)
266 {
267  if (n<k)
268  return(0);
269  gfan::Integer num = 1;
270  gfan::Integer den = 1;
271  for (int i=1; i<=k; i++)
272  den = den*i;
273  for (int j=n-k+1; j<=n; j++)
274  num = num*j;
275  gfan::Integer bin = num/den;
276  return(bin.toInt());
277 }
278 
279 
280 intvec* intToAface(unsigned int v0, int n, int k)
281 {
282  intvec* v = new intvec(k);
283  int j = 0;
284  for (int i=0; i<n; i++)
285  {
286  if (v0 & (1<<i))
287  (*v)[j++] = i+1;
288  }
289  return v;
290 }
291 
292 
294 {
295  leftv u = args;
296  if ((u != NULL) && (u->Typ() == INT_CMD))
297  {
298  leftv v = u->next;
299  if ((v != NULL) && (v->Typ() == INT_CMD))
300  {
301  int n = (int)(long) u->Data();
302  int k = (int)(long) v->Data();
303  unsigned int v = 0;
304  for (int i=0; i<k; i++)
305  v |= 1<<i; // sets the first k bits of v as 1
306 
308  int count = (int) binomial(n,k); L->Init(count);
309  unsigned int t;
310  while (!(v & (1<<n)))
311  {
312  L->m[--count].rtyp = INTVEC_CMD;
313  L->m[count].data = (void*) intToAface(v,n,k);
314 
315  // t gets v's least significant 0 bits set to 1
316  t = v | (v - 1);
317  // Next set to 1 the most significant bit to change,
318  // set to 0 the least significant ones, and add the necessary 1 bits.
319  v = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));
320  }
321  res->rtyp = LIST_CMD;
322  res->data = (void*) L;
323  return FALSE;
324  }
325  }
326  WerrorS("listOfAfacesToCheck: unexpected parameter");
327  return TRUE;
328 }
329 
330 
332 {
333  leftv u = args;
334  if ((u != NULL) && (u->Typ() == INTVEC_CMD))
335  {
336  leftv v = u->next;
337  if ((v != NULL) && (v->Typ() == INT_CMD))
338  {
339  leftv w = v->next;
340  if ((w != NULL) && (w->Typ() == INT_CMD))
341  {
342  intvec* aface = (intvec*) u->Data();
343  int ambientDimension = (int)(long) v->Data();
344  int dimension = (int)(long) w->Data();
345 
346  unsigned int af = 0;
347  for (int i=0; i<aface->length(); i++)
348  af |= 1<<((*aface)[i]-1);
349 
350  unsigned int t = af | (af - 1);
351  af = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(af) + 1));
352 
353  if (af & (1<<ambientDimension))
354  {
355  res->rtyp = INTVEC_CMD;
356  res->data = (void*) new intvec(1);
357  return FALSE;
358  }
359 
360  res->rtyp = INTVEC_CMD;
361  res->data = (void*) intToAface(af,ambientDimension,dimension);
362  return FALSE;
363  }
364  }
365  }
366  WerrorS("nextAfaceToCheck: unexpected parameter");
367  return TRUE;
368 }
369 
370 
372 {
373  p->iiAddCproc("","refineCones",FALSE,refineCones);
374  p->iiAddCproc("","listOfAfacesToCheck",FALSE,listOfAfacesToCheck);
375  p->iiAddCproc("","nextAfaceToCheck",FALSE,nextAfaceToCheck);
376 }
377 
378 #endif
int status int void size_t count
Definition: si_signals.h:59
bigintmat * transpose()
Definition: bigintmat.cc:38
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
static int binomial(int n, int k)
Definition: gitfan.cc:265
sleftv * m
Definition: lists.h:45
Class used for (list of) interpreter objects.
Definition: subexpr.h:82
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
void mu(int **points, int sizePoints)
Definition: tok.h:95
Definition: lists.h:22
CanonicalForm num(const CanonicalForm &f)
int check
Definition: libparse.cc:1104
#define FALSE
Definition: auxiliary.h:94
return P p
Definition: myNF.cc:203
Matrices of numbers.
Definition: bigintmat.h:51
Definition: gitfan.cc:28
gfan::ZVector facetNormal
Definition: gitfan.h:21
#define TRUE
Definition: auxiliary.h:98
gfan::ZCone eta
Definition: gitfan.h:19
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
int fanID
Definition: bbfan.cc:20
BOOLEAN listOfAfacesToCheck(leftv res, leftv args)
Definition: gitfan.cc:293
int Typ()
Definition: subexpr.cc:995
gfan::ZVector randomPoint(const gfan::ZCone *zc)
Definition: bbcone.cc:1056
void * data
Definition: subexpr.h:88
poly res
Definition: myNF.cc:322
void mergeFacets(facets &F, const facets &newFacets)
Definition: gitfan.cc:80
BOOLEAN inequalities(leftv res, leftv args)
Definition: bbcone.cc:561
const ring r
Definition: syzextra.cc:208
Definition: intvec.h:14
int j
Definition: myNF.cc:70
#define assume(x)
Definition: mod2.h:394
gfan::ZMatrix * bigintmatToZMatrix(const bigintmat &bim)
std::set< facet, facet_compare > facets
Definition: gitfan.h:50
BOOLEAN nextAfaceToCheck(leftv res, leftv args)
Definition: gitfan.cc:331
FILE * f
Definition: checklibs.c:9
int i
Definition: cfEzgcd.cc:123
Variable next() const
Definition: factory.h:135
leftv next
Definition: subexpr.h:86
INLINE_THIS void Init(int l=0)
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
int nr
Definition: lists.h:43
BOOLEAN dimension(leftv res, leftv args)
Definition: bbcone.cc:758
#define NULL
Definition: omList.c:10
gfan::ZVector interiorPoint
Definition: gitfan.h:20
slists * lists
Definition: mpr_numeric.h:146
int length() const
Definition: intvec.h:86
CanonicalForm den(const CanonicalForm &f)
static gfan::ZCone subcone(const lists &cones, const gfan::ZVector &point)
Definition: gitfan.cc:94
BOOLEAN equations(leftv res, leftv args)
Definition: bbcone.cc:578
const CanonicalForm & w
Definition: facAbsFact.cc:55
void lambda(int **points, int sizePoints)
int rtyp
Definition: subexpr.h:91
BOOLEAN ambientDimension(leftv res, leftv args)
Definition: bbcone.cc:724
void * Data()
Definition: subexpr.cc:1137
Definition: tok.h:117
void gitfan_setup(SModulFunctions *p)
Definition: gitfan.cc:371
omBin slists_bin
Definition: lists.cc:23
intvec * intToAface(unsigned int v0, int n, int k)
Definition: gitfan.cc:280
int BOOLEAN
Definition: auxiliary.h:85
BOOLEAN refineCones(leftv res, leftv args)
Definition: gitfan.cc:165
static gitfan::facets interiorFacets(const gfan::ZCone &zc, const gfan::ZCone &bound)
Definition: gitfan.cc:107