?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
rename cdd2polylib.pl to latte2polylib.pl
[barvinok.git] / util.c
blob57eecafd8c7c0cf9dd8fe9aa3184eabfe235ca81
1 #include <polylib/polylibgmp.h>
2 #include <stdlib.h>
3 #include <assert.h>
4 #include "config.h"
5 #include "version.h"
6
7 #ifndef HAVE_ENUMERATE4
8 #define Polyhedron_Enumerate(a,b,c,d) Polyhedron_Enumerate(a,b,c)
9 #endif
10
11 #define ALLOC(type) (type*)malloc(sizeof(type))
12 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
13
14 #ifdef __GNUC__
15 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
16 #else
17 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
18 #endif
19
20 #ifndef HAVE_ENUMERATION_FREE
21 #define Enumeration_Free(en) /* just leak some memory */
22 #endif
23
24 void manual_count(Polyhedron *P, Value* result)
25 {
26 Polyhedron *U = Universe_Polyhedron(0);
27 Enumeration *en = Polyhedron_Enumerate(P,U,1024,NULL);
28 Value *v = compute_poly(en,NULL);
29 value_assign(*result, *v);
30 value_clear(*v);
31 free(v);
32 Enumeration_Free(en);
33 Polyhedron_Free(U);
34 }
35
36 #ifndef HAVE_ENUMERATION_FREE
37 #undef Enumeration_Free
38 #endif
39
40 #include <barvinok/evalue.h>
41 #include <barvinok/util.h>
42 #include <barvinok/barvinok.h>
43
44 /* Return random value between 0 and max-1 inclusive
45 */
46 int random_int(int max) {
47 return (int) (((double)(max))*rand()/(RAND_MAX+1.0));
48 }
49
50 /* Inplace polarization
51 */
52 void Polyhedron_Polarize(Polyhedron *P)
53 {
54 unsigned NbRows = P->NbConstraints + P->NbRays;
55 int i;
56 Value **q;
57
58 q = (Value **)malloc(NbRows * sizeof(Value *));
59 assert(q);
60 for (i = 0; i < P->NbRays; ++i)
61 q[i] = P->Ray[i];
62 for (; i < NbRows; ++i)
63 q[i] = P->Constraint[i-P->NbRays];
64 P->NbConstraints = NbRows - P->NbConstraints;
65 P->NbRays = NbRows - P->NbRays;
66 free(P->Constraint);
67 P->Constraint = q;
68 P->Ray = q + P->NbConstraints;
69 }
70
71 /*
72 * Rather general polar
73 * We can optimize it significantly if we assume that
74 * P includes zero
75 *
76 * Also, we calculate the polar as defined in Schrijver
77 * The opposite should probably work as well and would
78 * eliminate the need for multiplying by -1
79 */
80 Polyhedron* Polyhedron_Polar(Polyhedron *P, unsigned NbMaxRays)
81 {
82 int i;
83 Value mone;
84 unsigned dim = P->Dimension + 2;
85 Matrix *M = Matrix_Alloc(P->NbRays, dim);
86
87 assert(M);
88 value_init(mone);
89 value_set_si(mone, -1);
90 for (i = 0; i < P->NbRays; ++i) {
91 Vector_Scale(P->Ray[i], M->p[i], mone, dim);
92 value_multiply(M->p[i][0], M->p[i][0], mone);
93 value_multiply(M->p[i][dim-1], M->p[i][dim-1], mone);
94 }
95 P = Constraints2Polyhedron(M, NbMaxRays);
96 assert(P);
97 Matrix_Free(M);
98 value_clear(mone);
99 return P;
100 }
101
102 /*
103 * Returns the supporting cone of P at the vertex with index v
104 */
105 Polyhedron* supporting_cone(Polyhedron *P, int v)
106 {
107 Matrix *M;
108 Value tmp;
109 int i, n, j;
110 unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
111 unsigned dim = P->Dimension + 2;
112
113 assert(v >=0 && v < P->NbRays);
114 assert(value_pos_p(P->Ray[v][dim-1]));
115 assert(supporting);
116
117 value_init(tmp);
118 for (i = 0, n = 0; i < P->NbConstraints; ++i) {
119 Inner_Product(P->Constraint[i] + 1, P->Ray[v] + 1, dim - 1, &tmp);
120 if ((supporting[i] = value_zero_p(tmp)))
121 ++n;
122 }
123 assert(n >= dim - 2);
124 value_clear(tmp);
125 M = Matrix_Alloc(n, dim);
126 assert(M);
127 for (i = 0, j = 0; i < P->NbConstraints; ++i)
128 if (supporting[i]) {
129 value_set_si(M->p[j][dim-1], 0);
130 Vector_Copy(P->Constraint[i], M->p[j++], dim-1);
131 }
132 free(supporting);
133 P = Constraints2Polyhedron(M, P->NbRays+1);
134 assert(P);
135 Matrix_Free(M);
136 return P;
137 }
138
139 void value_lcm(Value i, Value j, Value* lcm)
140 {
141 Value aux;
142 value_init(aux);
143 value_multiply(aux,i,j);
144 Gcd(i,j,lcm);
145 value_division(*lcm,aux,*lcm);
146 value_clear(aux);
147 }
148
149 Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v)
150 {
151 Matrix *M;
152 Value lcm, tmp, tmp2;
153 unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
154 unsigned dim = P->Dimension + 2;
155 unsigned nparam = v->Vertex->NbColumns - 2;
156 unsigned nvar = dim - nparam - 2;
157 int i, n, j;
158 Vector *row;
159
160 assert(supporting);
161 row = Vector_Alloc(nparam+1);
162 assert(row);
163 value_init(lcm);
164 value_init(tmp);
165 value_init(tmp2);
166 value_set_si(lcm, 1);
167 for (i = 0, n = 0; i < P->NbConstraints; ++i) {
168 Vector_Set(row->p, 0, nparam+1);
169 for (j = 0 ; j < nvar; ++j) {
170 value_set_si(tmp, 1);
171 value_assign(tmp2, P->Constraint[i][j+1]);
172 if (value_ne(lcm, v->Vertex->p[j][nparam+1])) {
173 value_assign(tmp, lcm);
174 value_lcm(lcm, v->Vertex->p[j][nparam+1], &lcm);
175 value_division(tmp, lcm, tmp);
176 value_multiply(tmp2, tmp2, lcm);
177 value_division(tmp2, tmp2, v->Vertex->p[j][nparam+1]);
178 }
179 Vector_Combine(row->p, v->Vertex->p[j], row->p,
180 tmp, tmp2, nparam+1);
181 }
182 value_set_si(tmp, 1);
183 Vector_Combine(row->p, P->Constraint[i]+1+nvar, row->p, tmp, lcm, nparam+1);
184 for (j = 0; j < nparam+1; ++j)
185 if (value_notzero_p(row->p[j]))
186 break;
187 if ((supporting[i] = (j == nparam + 1)))
188 ++n;
189 }
190 assert(n >= nvar);
191 value_clear(tmp);
192 value_clear(tmp2);
193 value_clear(lcm);
194 Vector_Free(row);
195 M = Matrix_Alloc(n, nvar+2);
196 assert(M);
197 for (i = 0, j = 0; i < P->NbConstraints; ++i)
198 if (supporting[i]) {
199 value_set_si(M->p[j][nvar+1], 0);
200 Vector_Copy(P->Constraint[i], M->p[j++], nvar+1);
201 }
202 free(supporting);
203 P = Constraints2Polyhedron(M, P->NbRays+1);
204 assert(P);
205 Matrix_Free(M);
206 return P;
207 }
208
209 Polyhedron* triangulate_cone(Polyhedron *P, unsigned NbMaxCons)
210 {
211 const static int MAX_TRY=10;
212 int i, j, r, n, t;
213 Value tmp;
214 unsigned dim = P->Dimension;
215 Matrix *M = Matrix_Alloc(P->NbRays+1, dim+3);
216 Matrix *M2, *M3;
217 Polyhedron *L, *R, *T;
218 assert(P->NbEq == 0);
219
220 R = NULL;
221 value_init(tmp);
222
223 Vector_Set(M->p[0]+1, 0, dim+1);
224 value_set_si(M->p[0][0], 1);
225 value_set_si(M->p[0][dim+2], 1);
226 Vector_Set(M->p[P->NbRays]+1, 0, dim+2);
227 value_set_si(M->p[P->NbRays][0], 1);
228 value_set_si(M->p[P->NbRays][dim+1], 1);
229
230 /* Delaunay triangulation */
231 for (i = 0, r = 1; i < P->NbRays; ++i) {
232 if (value_notzero_p(P->Ray[i][dim+1]))
233 continue;
234 Vector_Copy(P->Ray[i], M->p[r], dim+1);
235 Inner_Product(M->p[r]+1, M->p[r]+1, dim, &tmp);
236 value_assign(M->p[r][dim+1], tmp);
237 value_set_si(M->p[r][dim+2], 0);
238 ++r;
239 }
240
241 M3 = Matrix_Copy(M);
242 L = Rays2Polyhedron(M3, NbMaxCons);
243 Matrix_Free(M3);
244
245 M2 = Matrix_Alloc(dim+1, dim+2);
246
247 t = 1;
248 if (0) {
249 try_again:
250 /* Usually R should still be 0 */
251 Domain_Free(R);
252 Polyhedron_Free(L);
253 for (r = 1; r < P->NbRays; ++r) {
254 value_set_si(M->p[r][dim+1], random_int((t+1)*dim)+1);
255 }
256 M3 = Matrix_Copy(M);
257 L = Rays2Polyhedron(M3, NbMaxCons);
258 Matrix_Free(M3);
259 ++t;
260 }
261 assert(t <= MAX_TRY);
262
263 R = NULL;
264 n = 0;
265
266 for (i = 0; i < L->NbConstraints; ++i) {
267 /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
268 if (value_negz_p(L->Constraint[i][dim+1]))
269 continue;
270 if (value_notzero_p(L->Constraint[i][dim+2]))
271 continue;
272 for (j = 1, r = 1; j < M->NbRows; ++j) {
273 Inner_Product(M->p[j]+1, L->Constraint[i]+1, dim+1, &tmp);
274 if (value_notzero_p(tmp))
275 continue;
276 if (r > dim)
277 goto try_again;
278 Vector_Copy(M->p[j]+1, M2->p[r]+1, dim);
279 value_set_si(M2->p[r][0], 1);
280 value_set_si(M2->p[r][dim+1], 0);
281 ++r;
282 }
283 assert(r == dim+1);
284 Vector_Set(M2->p[0]+1, 0, dim);
285 value_set_si(M2->p[0][0], 1);
286 value_set_si(M2->p[0][dim+1], 1);
287 T = Rays2Polyhedron(M2, P->NbConstraints+1);
288 T->next = R;
289 R = T;
290 ++n;
291 }
292 Matrix_Free(M2);
293
294 Polyhedron_Free(L);
295 value_clear(tmp);
296 Matrix_Free(M);
297
298 return R;
299 }
300
301 void check_triangulization(Polyhedron *P, Polyhedron *T)
302 {
303 Polyhedron *C, *D, *E, *F, *G, *U;
304 for (C = T; C; C = C->next) {
305 if (C == T)
306 U = C;
307 else
308 U = DomainConvex(DomainUnion(U, C, 100), 100);
309 for (D = C->next; D; D = D->next) {
310 F = C->next;
311 G = D->next;
312 C->next = NULL;
313 D->next = NULL;
314 E = DomainIntersection(C, D, 600);
315 assert(E->NbRays == 0 || E->NbEq >= 1);
316 Polyhedron_Free(E);
317 C->next = F;
318 D->next = G;
319 }
320 }
321 assert(PolyhedronIncludes(U, P));
322 assert(PolyhedronIncludes(P, U));
323 }
324
325 /* Computes x, y and g such that g = gcd(a,b) and a*x+b*y = g */
326 static void Euclid(Value a, Value b, Value *x, Value *y, Value *g)
327 {
328 Value c, d, e, f, tmp;
329
330 value_init(c);
331 value_init(d);
332 value_init(e);
333 value_init(f);
334 value_init(tmp);
335 value_absolute(c, a);
336 value_absolute(d, b);
337 value_set_si(e, 1);
338 value_set_si(f, 0);
339 while(value_pos_p(d)) {
340 value_division(tmp, c, d);
341 value_multiply(tmp, tmp, f);
342 value_subtract(e, e, tmp);
343 value_division(tmp, c, d);
344 value_multiply(tmp, tmp, d);
345 value_subtract(c, c, tmp);
346 value_swap(c, d);
347 value_swap(e, f);
348 }
349 value_assign(*g, c);
350 if (value_zero_p(a))
351 value_set_si(*x, 0);
352 else if (value_pos_p(a))
353 value_assign(*x, e);
354 else value_oppose(*x, e);
355 if (value_zero_p(b))
356 value_set_si(*y, 0);
357 else {
358 value_multiply(tmp, a, *x);
359 value_subtract(tmp, c, tmp);
360 value_division(*y, tmp, b);
361 }
362 value_clear(c);
363 value_clear(d);
364 value_clear(e);
365 value_clear(f);
366 value_clear(tmp);
367 }
368
369 Matrix * unimodular_complete(Vector *row)
370 {
371 Value g, b, c, old, tmp;
372 Matrix *m;
373 unsigned i, j;
374
375 value_init(b);
376 value_init(c);
377 value_init(g);
378 value_init(old);
379 value_init(tmp);
380 m = Matrix_Alloc(row->Size, row->Size);
381 for (j = 0; j < row->Size; ++j) {
382 value_assign(m->p[0][j], row->p[j]);
383 }
384 value_assign(g, row->p[0]);
385 for (i = 1; value_zero_p(g) && i < row->Size; ++i) {
386 for (j = 0; j < row->Size; ++j) {
387 if (j == i-1)
388 value_set_si(m->p[i][j], 1);
389 else
390 value_set_si(m->p[i][j], 0);
391 }
392 value_assign(g, row->p[i]);
393 }
394 for (; i < row->Size; ++i) {
395 value_assign(old, g);
396 Euclid(old, row->p[i], &c, &b, &g);
397 value_oppose(b, b);
398 for (j = 0; j < row->Size; ++j) {
399 if (j < i) {
400 value_multiply(tmp, row->p[j], b);
401 value_division(m->p[i][j], tmp, old);
402 } else if (j == i)
403 value_assign(m->p[i][j], c);
404 else
405 value_set_si(m->p[i][j], 0);
406 }
407 }
408 value_clear(b);
409 value_clear(c);
410 value_clear(g);
411 value_clear(old);
412 value_clear(tmp);
413 return m;
414 }
415
416 /*
417 * Returns a full-dimensional polyhedron with the same number
418 * of integer points as P
419 */
420 Polyhedron *remove_equalities(Polyhedron *P)
421 {
422 Value g;
423 Vector *v;
424 Polyhedron *p = Polyhedron_Copy(P), *q;
425 unsigned dim = p->Dimension;
426 Matrix *m1, *m2;
427 int i;
428
429 value_init(g);
430 while (p->NbEq > 0) {
431 assert(dim > 0);
432 Vector_Gcd(p->Constraint[0]+1, dim+1, &g);
433 Vector_AntiScale(p->Constraint[0]+1, p->Constraint[0]+1, g, dim+1);
434 Vector_Gcd(p->Constraint[0]+1, dim, &g);
435 if (value_notone_p(g) && value_notmone_p(g)) {
436 Polyhedron_Free(p);
437 p = Empty_Polyhedron(0);
438 break;
439 }
440 v = Vector_Alloc(dim);
441 Vector_Copy(p->Constraint[0]+1, v->p, dim);
442 m1 = unimodular_complete(v);
443 m2 = Matrix_Alloc(dim, dim+1);
444 for (i = 0; i < dim-1 ; ++i) {
445 Vector_Copy(m1->p[i+1], m2->p[i], dim);
446 value_set_si(m2->p[i][dim], 0);
447 }
448 Vector_Set(m2->p[dim-1], 0, dim);
449 value_set_si(m2->p[dim-1][dim], 1);
450 q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
451 Vector_Free(v);
452 Matrix_Free(m1);
453 Matrix_Free(m2);
454 Polyhedron_Free(p);
455 p = q;
456 --dim;
457 }
458 value_clear(g);
459 return p;
460 }
461
462 /*
463 * Returns a full-dimensional polyhedron with the same number
464 * of integer points as P
465 * nvar specifies the number of variables
466 * The remaining dimensions are assumed to be parameters
467 * Destroys P
468 * factor is NbEq x (nparam+2) matrix, containing stride constraints
469 * on the parameters; column nparam is the constant;
470 * column nparam+1 is the stride
471 *
472 * if factor is NULL, only remove equalities that don't affect
473 * the number of points
474 */
475 Polyhedron *remove_equalities_p(Polyhedron *P, unsigned nvar, Matrix **factor)
476 {
477 Value g;
478 Vector *v;
479 Polyhedron *p = P, *q;
480 unsigned dim = p->Dimension;
481 Matrix *m1, *m2, *f;
482 int i, j, skip;
483
484 value_init(g);
485 if (factor) {
486 f = Matrix_Alloc(p->NbEq, dim-nvar+2);
487 *factor = f;
488 }
489 j = 0;
490 skip = 0;
491 while (nvar > 0 && p->NbEq - skip > 0) {
492 assert(dim > 0);
493
494 while (value_zero_p(p->Constraint[skip][0]) &&
495 First_Non_Zero(p->Constraint[skip]+1, nvar) == -1)
496 ++skip;
497 if (p->NbEq == skip)
498 break;
499
500 Vector_Gcd(p->Constraint[skip]+1, dim+1, &g);
501 Vector_AntiScale(p->Constraint[skip]+1, p->Constraint[skip]+1, g, dim+1);
502 Vector_Gcd(p->Constraint[skip]+1, nvar, &g);
503 if (!factor && value_notone_p(g) && value_notmone_p(g)) {
504 ++skip;
505 continue;
506 }
507 if (factor) {
508 Vector_Copy(p->Constraint[skip]+1+nvar, f->p[j], dim-nvar+1);
509 value_assign(f->p[j][dim-nvar+1], g);
510 }
511 v = Vector_Alloc(dim);
512 Vector_AntiScale(p->Constraint[skip]+1, v->p, g, nvar);
513 Vector_Set(v->p+nvar, 0, dim-nvar);
514 m1 = unimodular_complete(v);
515 m2 = Matrix_Alloc(dim, dim+1);
516 for (i = 0; i < dim-1 ; ++i) {
517 Vector_Copy(m1->p[i+1], m2->p[i], dim);
518 value_set_si(m2->p[i][dim], 0);
519 }
520 Vector_Set(m2->p[dim-1], 0, dim);
521 value_set_si(m2->p[dim-1][dim], 1);
522 q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
523 Vector_Free(v);
524 Matrix_Free(m1);
525 Matrix_Free(m2);
526 Polyhedron_Free(p);
527 p = q;
528 --dim;
529 --nvar;
530 ++j;
531 }
532 value_clear(g);
533 return p;
534 }
535
536 void Line_Length(Polyhedron *P, Value *len)
537 {
538 Value tmp, pos, neg;
539 int p = 0, n = 0;
540 int i;
541
542 assert(P->Dimension == 1);
543
544 value_init(tmp);
545 value_init(pos);
546 value_init(neg);
547
548 for (i = 0; i < P->NbConstraints; ++i) {
549 value_oppose(tmp, P->Constraint[i][2]);
550 if (value_pos_p(P->Constraint[i][1])) {
551 mpz_cdiv_q(tmp, tmp, P->Constraint[i][1]);
552 if (!p || value_gt(tmp, pos))
553 value_assign(pos, tmp);
554 p = 1;
555 } else {
556 mpz_fdiv_q(tmp, tmp, P->Constraint[i][1]);
557 if (!n || value_lt(tmp, neg))
558 value_assign(neg, tmp);
559 n = 1;
560 }
561 if (n && p) {
562 value_subtract(tmp, neg, pos);
563 value_increment(*len, tmp);
564 } else
565 value_set_si(*len, -1);
566 }
567
568 value_clear(tmp);
569 value_clear(pos);
570 value_clear(neg);
571 }
572
573 /*
574 * Factors the polyhedron P into polyhedra Q_i such that
575 * the number of integer points in P is equal to the product
576 * of the number of integer points in the individual Q_i
577 *
578 * If no factors can be found, NULL is returned.
579 * Otherwise, a linked list of the factors is returned.
580 *
581 * The algorithm works by first computing the Hermite normal form
582 * and then grouping columns linked by one or more constraints together,
583 * where a constraints "links" two or more columns if the constraint
584 * has nonzero coefficients in the columns.
585 */
586 Polyhedron* Polyhedron_Factor(Polyhedron *P, unsigned nparam,
587 unsigned NbMaxRays)
588 {
589 int i, j, k;
590 Matrix *M, *H, *Q, *U;
591 int *pos; /* for each column: row position of pivot */
592 int *group; /* group to which a column belongs */
593 int *cnt; /* number of columns in the group */
594 int *rowgroup; /* group to which a constraint belongs */
595 int nvar = P->Dimension - nparam;
596 Polyhedron *F = NULL;
597
598 if (nvar <= 1)
599 return NULL;
600
601 NALLOC(pos, nvar);
602 NALLOC(group, nvar);
603 NALLOC(cnt, nvar);
604 NALLOC(rowgroup, P->NbConstraints);
605
606 M = Matrix_Alloc(P->NbConstraints, nvar);
607 for (i = 0; i < P->NbConstraints; ++i)
608 Vector_Copy(P->Constraint[i]+1, M->p[i], nvar);
609 left_hermite(M, &H, &Q, &U);
610 Matrix_Free(M);
611 Matrix_Free(Q);
612 Matrix_Free(U);
613
614 for (i = 0; i < P->NbConstraints; ++i)
615 rowgroup[i] = -1;
616 for (i = 0, j = 0; i < H->NbColumns; ++i) {
617 for ( ; j < H->NbRows; ++j)
618 if (value_notzero_p(H->p[j][i]))
619 break;
620 assert (j < H->NbRows);
621 pos[i] = j;
622 }
623 for (i = 0; i < nvar; ++i) {
624 group[i] = i;
625 cnt[i] = 1;
626 }
627 for (i = 0; i < H->NbColumns && cnt[0] < nvar; ++i) {
628 if (rowgroup[pos[i]] == -1)
629 rowgroup[pos[i]] = i;
630 for (j = pos[i]+1; j < H->NbRows; ++j) {
631 if (value_zero_p(H->p[j][i]))
632 continue;
633 if (rowgroup[j] != -1)
634 continue;
635 rowgroup[j] = group[i];
636 for (k = i+1; k < H->NbColumns && j >= pos[k]; ++k) {
637 int g = group[k];
638 while (cnt[g] == 0)
639 g = group[g];
640 group[k] = g;
641 if (group[k] != group[i] && value_notzero_p(H->p[j][k])) {
642 assert(cnt[group[k]] != 0);
643 assert(cnt[group[i]] != 0);
644 if (group[i] < group[k]) {
645 cnt[group[i]] += cnt[group[k]];
646 cnt[group[k]] = 0;
647 group[k] = group[i];
648 } else {
649 cnt[group[k]] += cnt[group[i]];
650 cnt[group[i]] = 0;
651 group[i] = group[k];
652 }
653 }
654 }
655 }
656 }
657
658 if (cnt[0] != nvar) {
659 /* Extract out pure context constraints separately */
660 Polyhedron **next = &F;
661 for (i = nparam ? -1 : 0; i < nvar; ++i) {
662 int d;
663
664 if (i == -1) {
665 for (j = 0, k = 0; j < P->NbConstraints; ++j)
666 if (rowgroup[j] == -1) {
667 if (First_Non_Zero(P->Constraint[j]+1+nvar,
668 nparam) == -1)
669 rowgroup[j] = -2;
670 else
671 ++k;
672 }
673 if (k == 0)
674 continue;
675 d = 0;
676 } else {
677 if (cnt[i] == 0)
678 continue;
679 d = cnt[i];
680 for (j = 0, k = 0; j < P->NbConstraints; ++j)
681 if (rowgroup[j] >= 0 && group[rowgroup[j]] == i) {
682 rowgroup[j] = i;
683 ++k;
684 }
685 }
686
687 M = Matrix_Alloc(k, d+nparam+2);
688 for (j = 0, k = 0; j < P->NbConstraints; ++j) {
689 int l, m;
690 if (rowgroup[j] != i)
691 continue;
692 value_assign(M->p[k][0], P->Constraint[j][0]);
693 for (l = 0, m = 0; m < d; ++l) {
694 if (group[l] != i)
695 continue;
696 value_assign(M->p[k][1+m++], H->p[j][l]);
697 }
698 Vector_Copy(P->Constraint[j]+1+nvar, M->p[k]+1+m, nparam+1);
699 ++k;
700 }
701 *next = Constraints2Polyhedron(M, NbMaxRays);
702 next = &(*next)->next;
703 Matrix_Free(M);
704 }
705 }
706 Matrix_Free(H);
707 free(pos);
708 free(group);
709 free(cnt);
710 free(rowgroup);
711 return F;
712 }
713
714 /*
715 * Project on final dim dimensions
716 */
717 Polyhedron* Polyhedron_Project(Polyhedron *P, int dim)
718 {
719 int i;
720 int remove = P->Dimension - dim;
721 Matrix *T;
722 Polyhedron *I;
723
724 if (P->Dimension == dim)
725 return Polyhedron_Copy(P);
726
727 T = Matrix_Alloc(dim+1, P->Dimension+1);
728 for (i = 0; i < dim+1; ++i)
729 value_set_si(T->p[i][i+remove], 1);
730 I = Polyhedron_Image(P, T, P->NbConstraints);
731 Matrix_Free(T);
732 return I;
733 }
734
735 /* Constructs a new constraint that ensures that
736 * the first constraint is (strictly) smaller than
737 * the second.
738 */
739 static void smaller_constraint(Value *a, Value *b, Value *c, int pos, int shift,
740 int len, int strict, Value *tmp)
741 {
742 value_oppose(*tmp, b[pos+1]);
743 value_set_si(c[0], 1);
744 Vector_Combine(a+1+shift, b+1+shift, c+1, *tmp, a[pos+1], len-shift-1);
745 if (strict)
746 value_decrement(c[len-shift-1], c[len-shift-1]);
747 ConstraintSimplify(c, c, len-shift, tmp);
748 }
749
750 struct section { Polyhedron * D; evalue E; };
751
752 evalue * ParamLine_Length_mod(Polyhedron *P, Polyhedron *C, int MaxRays)
753 {
754 unsigned dim = P->Dimension;
755 unsigned nvar = dim - C->Dimension;
756 int *pos;
757 int i, j, p, n, z;
758 struct section *s;
759 Matrix *M, *M2;
760 int nd = 0;
761 int k, l, k2, l2, q;
762 evalue *L, *U;
763 evalue *F;
764 Value g;
765 Polyhedron *T;
766 evalue mone;
767
768 assert(nvar == 1);
769
770 NALLOC(pos, P->NbConstraints);
771 value_init(g);
772 value_init(mone.d);
773 evalue_set_si(&mone, -1, 1);
774
775 for (i = 0, z = 0; i < P->NbConstraints; ++i)
776 if (value_zero_p(P->Constraint[i][1]))
777 ++z;
778 /* put those with positive coefficients first; number: p */
779 for (i = 0, p = 0, n = P->NbConstraints-z-1; i < P->NbConstraints; ++i)
780 if (value_pos_p(P->Constraint[i][1]))
781 pos[p++] = i;
782 else if (value_neg_p(P->Constraint[i][1]))
783 pos[n--] = i;
784 n = P->NbConstraints-z-p;
785 assert (p >= 1 && n >= 1);
786 s = (struct section *) malloc(p * n * sizeof(struct section));
787 M = Matrix_Alloc((p-1) + (n-1), dim-nvar+2);
788 for (k = 0; k < p; ++k) {
789 for (k2 = 0; k2 < p; ++k2) {
790 if (k2 == k)
791 continue;
792 q = k2 - (k2 > k);
793 smaller_constraint(
794 P->Constraint[pos[k]],
795 P->Constraint[pos[k2]],
796 M->p[q], 0, nvar, dim+2, k2 > k, &g);
797 }
798 for (l = p; l < p+n; ++l) {
799 for (l2 = p; l2 < p+n; ++l2) {
800 if (l2 == l)
801 continue;
802 q = l2-1 - (l2 > l);
803 smaller_constraint(
804 P->Constraint[pos[l2]],
805 P->Constraint[pos[l]],
806 M->p[q], 0, nvar, dim+2, l2 > l, &g);
807 }
808 M2 = Matrix_Copy(M);
809 T = Constraints2Polyhedron(M2, P->NbRays);
810 Matrix_Free(M2);
811 s[nd].D = DomainIntersection(T, C, MaxRays);
812 Domain_Free(T);
813 POL_ENSURE_VERTICES(s[nd].D);
814 if (emptyQ(s[nd].D)) {
815 Polyhedron_Free(s[nd].D);
816 continue;
817 }
818 L = bv_ceil3(P->Constraint[pos[k]]+1+nvar,
819 dim-nvar+1,
820 P->Constraint[pos[k]][0+1], s[nd].D);
821 U = bv_ceil3(P->Constraint[pos[l]]+1+nvar,
822 dim-nvar+1,
823 P->Constraint[pos[l]][0+1], s[nd].D);
824 eadd(L, U);
825 eadd(&mone, U);
826 emul(&mone, U);
827 s[nd].E = *U;
828 free_evalue_refs(L);
829 free(L);
830 free(U);
831 ++nd;
832 }
833 }
834
835 Matrix_Free(M);
836
837 F = ALLOC(evalue);
838 value_init(F->d);
839 value_set_si(F->d, 0);
840 F->x.p = new_enode(partition, 2*nd, dim-nvar);
841 for (k = 0; k < nd; ++k) {
842 EVALUE_SET_DOMAIN(F->x.p->arr[2*k], s[k].D);
843 value_clear(F->x.p->arr[2*k+1].d);
844 F->x.p->arr[2*k+1] = s[k].E;
845 }
846 free(s);
847
848 free_evalue_refs(&mone);
849 value_clear(g);
850 free(pos);
851
852 return F;
853 }
854
855 #ifdef USE_MODULO
856 evalue* ParamLine_Length(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
857 {
858 return ParamLine_Length_mod(P, C, MaxRays);
859 }
860 #else
861 evalue* ParamLine_Length(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
862 {
863 evalue* tmp;
864 tmp = ParamLine_Length_mod(P, C, MaxRays);
865 evalue_mod2table(tmp, C->Dimension);
866 reduce_evalue(tmp);
867 return tmp;
868 }
869 #endif
870
871 Bool isIdentity(Matrix *M)
872 {
873 unsigned i, j;
874 if (M->NbRows != M->NbColumns)
875 return False;
876
877 for (i = 0;i < M->NbRows; i ++)
878 for (j = 0; j < M->NbColumns; j ++)
879 if (i == j) {
880 if(value_notone_p(M->p[i][j]))
881 return False;
882 } else {
883 if(value_notzero_p(M->p[i][j]))
884 return False;
885 }
886 return True;
887 }
888
889 void Param_Polyhedron_Print(FILE* DST, Param_Polyhedron *PP, char **param_names)
890 {
891 Param_Domain *P;
892 Param_Vertices *V;
893
894 for(P=PP->D;P;P=P->next) {
895
896 /* prints current val. dom. */
897 printf( "---------------------------------------\n" );
898 printf( "Domain :\n");
899 Print_Domain( stdout, P->Domain, param_names );
900
901 /* scan the vertices */
902 printf( "Vertices :\n");
903 FORALL_PVertex_in_ParamPolyhedron(V,P,PP) {
904
905 /* prints each vertex */
906 Print_Vertex( stdout, V->Vertex, param_names );
907 printf( "\n" );
908 }
909 END_FORALL_PVertex_in_ParamPolyhedron;
910 }
911 }
912
913 void Enumeration_Print(FILE *Dst, Enumeration *en, char **params)
914 {
915 for (; en; en = en->next) {
916 Print_Domain(Dst, en->ValidityDomain, params);
917 print_evalue(Dst, &en->EP, params);
918 }
919 }
920
921 void Enumeration_Free(Enumeration *en)
922 {
923 Enumeration *ee;
924
925 while( en )
926 {
927 free_evalue_refs( &(en->EP) );
928 Domain_Free( en->ValidityDomain );
929 ee = en ->next;
930 free( en );
931 en = ee;
932 }
933 }
934
935 void Enumeration_mod2table(Enumeration *en, unsigned nparam)
936 {
937 for (; en; en = en->next) {
938 evalue_mod2table(&en->EP, nparam);
939 reduce_evalue(&en->EP);
940 }
941 }
942
943 size_t Enumeration_size(Enumeration *en)
944 {
945 size_t s = 0;
946
947 for (; en; en = en->next) {
948 s += domain_size(en->ValidityDomain);
949 s += evalue_size(&en->EP);
950 }
951 return s;
952 }
953
954 void Free_ParamNames(char **params, int m)
955 {
956 while (--m >= 0)
957 free(params[m]);
958 free(params);
959 }
960
961 int DomainIncludes(Polyhedron *Pol1, Polyhedron *Pol2)
962 {
963 Polyhedron *P2;
964 for ( ; Pol1; Pol1 = Pol1->next) {
965 for (P2 = Pol2; P2; P2 = P2->next)
966 if (!PolyhedronIncludes(Pol1, P2))
967 break;
968 if (!P2)
969 return 1;
970 }
971 return 0;
972 }
973
974 int line_minmax(Polyhedron *I, Value *min, Value *max)
975 {
976 int i;
977
978 if (I->NbEq >= 1) {
979 value_oppose(I->Constraint[0][2], I->Constraint[0][2]);
980 /* There should never be a remainder here */
981 if (value_pos_p(I->Constraint[0][1]))
982 mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
983 else
984 mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
985 value_assign(*max, *min);
986 } else for (i = 0; i < I->NbConstraints; ++i) {
987 if (value_zero_p(I->Constraint[i][1])) {
988 Polyhedron_Free(I);
989 return 0;
990 }
991
992 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
993 if (value_pos_p(I->Constraint[i][1]))
994 mpz_cdiv_q(*min, I->Constraint[i][2], I->Constraint[i][1]);
995 else
996 mpz_fdiv_q(*max, I->Constraint[i][2], I->Constraint[i][1]);
997 }
998 Polyhedron_Free(I);
999 return 1;
1000 }
1001
1002 /**
1003
1004 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1005 each imbriquation
1006
1007 @param pos index position of current loop index (1..hdim-1)
1008 @param P loop domain
1009 @param context context values for fixed indices
1010 @param exist number of existential variables
1011 @return the number of integer points in this
1012 polyhedron
1013
1014 */
1015 void count_points_e (int pos, Polyhedron *P, int exist, int nparam,
1016 Value *context, Value *res)
1017 {
1018 Value LB, UB, k, c;
1019
1020 if (emptyQ(P)) {
1021 value_set_si(*res, 0);
1022 return;
1023 }
1024
1025 value_init(LB); value_init(UB); value_init(k);
1026 value_set_si(LB,0);
1027 value_set_si(UB,0);
1028
1029 if (lower_upper_bounds(pos,P,context,&LB,&UB) !=0) {
1030 /* Problem if UB or LB is INFINITY */
1031 value_clear(LB); value_clear(UB); value_clear(k);
1032 if (pos > P->Dimension - nparam - exist)
1033 value_set_si(*res, 1);
1034 else
1035 value_set_si(*res, -1);
1036 return;
1037 }
1038
1039 #ifdef EDEBUG1
1040 if (!P->next) {
1041 int i;
1042 for (value_assign(k,LB); value_le(k,UB); value_increment(k,k)) {
1043 fprintf(stderr, "(");
1044 for (i=1; i<pos; i++) {
1045 value_print(stderr,P_VALUE_FMT,context[i]);
1046 fprintf(stderr,",");
1047 }
1048 value_print(stderr,P_VALUE_FMT,k);
1049 fprintf(stderr,")\n");
1050 }
1051 }
1052 #endif
1053
1054 value_set_si(context[pos],0);
1055 if (value_lt(UB,LB)) {
1056 value_clear(LB); value_clear(UB); value_clear(k);
1057 value_set_si(*res, 0);
1058 return;
1059 }
1060 if (!P->next) {
1061 if (exist)
1062 value_set_si(*res, 1);
1063 else {
1064 value_subtract(k,UB,LB);
1065 value_add_int(k,k,1);
1066 value_assign(*res, k);
1067 }
1068 value_clear(LB); value_clear(UB); value_clear(k);
1069 return;
1070 }
1071
1072 /*-----------------------------------------------------------------*/
1073 /* Optimization idea */
1074 /* If inner loops are not a function of k (the current index) */
1075 /* i.e. if P->Constraint[i][pos]==0 for all P following this and */
1076 /* for all i, */
1077 /* Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context) */
1078 /* (skip the for loop) */
1079 /*-----------------------------------------------------------------*/
1080
1081 value_init(c);
1082 value_set_si(*res, 0);
1083 for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
1084 /* Insert k in context */
1085 value_assign(context[pos],k);
1086 count_points_e(pos+1, P->next, exist, nparam, context, &c);
1087 if(value_notmone_p(c))
1088 value_addto(*res, *res, c);
1089 else {
1090 value_set_si(*res, -1);
1091 break;
1092 }
1093 if (pos > P->Dimension - nparam - exist &&
1094 value_pos_p(*res))
1095 break;
1096 }
1097 value_clear(c);
1098
1099 #ifdef EDEBUG11
1100 fprintf(stderr,"%d\n",CNT);
1101 #endif
1102
1103 /* Reset context */
1104 value_set_si(context[pos],0);
1105 value_clear(LB); value_clear(UB); value_clear(k);
1106 return;
1107 } /* count_points_e */
1108
1109 int DomainContains(Polyhedron *P, Value *list_args, int len,
1110 unsigned MaxRays, int set)
1111 {
1112 int i;
1113 Value m;
1114
1115 if (P->Dimension == len)
1116 return in_domain(P, list_args);
1117
1118 assert(set); // assume list_args is large enough
1119 assert((P->Dimension - len) % 2 == 0);
1120 value_init(m);
1121 for (i = 0; i < P->Dimension - len; i += 2) {
1122 int j, k;
1123 for (j = 0 ; j < P->NbEq; ++j)
1124 if (value_notzero_p(P->Constraint[j][1+len+i]))
1125 break;
1126 assert(j < P->NbEq);
1127 value_absolute(m, P->Constraint[j][1+len+i]);
1128 k = First_Non_Zero(P->Constraint[j]+1, len);
1129 assert(k != -1);
1130 assert(First_Non_Zero(P->Constraint[j]+1+k+1, len - k - 1) == -1);
1131 mpz_fdiv_q(list_args[len+i], list_args[k], m);
1132 mpz_fdiv_r(list_args[len+i+1], list_args[k], m);
1133 }
1134 value_clear(m);
1135
1136 return in_domain(P, list_args);
1137 }
1138
1139 Polyhedron *DomainConcat(Polyhedron *head, Polyhedron *tail)
1140 {
1141 Polyhedron *S;
1142 if (!head)
1143 return tail;
1144 for (S = head; S->next; S = S->next)
1145 ;
1146 S->next = tail;
1147 return head;
1148 }
1149
1150 #ifdef HAVE_LEXSMALLER
1151 #include <polylib/ranking.h>
1152
1153 evalue *barvinok_lexsmaller_ev(Polyhedron *P, Polyhedron *D, unsigned dim,
1154 Polyhedron *C, unsigned MaxRays)
1155 {
1156 evalue *ranking;
1157 Polyhedron *RC, *RD, *Q;
1158 unsigned nparam = dim + C->Dimension;
1159 unsigned exist;
1160 Polyhedron *CA;
1161
1162 RC = LexSmaller(P, D, dim, C, MaxRays);
1163 RD = RC->next;
1164 RC->next = NULL;
1165
1166 exist = RD->Dimension - nparam - dim;
1167 CA = align_context(RC, RD->Dimension, MaxRays);
1168 Q = DomainIntersection(RD, CA, MaxRays);
1169 Polyhedron_Free(CA);
1170 Domain_Free(RD);
1171 Polyhedron_Free(RC);
1172 RD = Q;
1173
1174 for (Q = RD; Q; Q = Q->next) {
1175 evalue *t;
1176 Polyhedron *next = Q->next;
1177 Q->next = 0;
1178
1179 t = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
1180
1181 if (Q == RD)
1182 ranking = t;
1183 else {
1184 eadd(t, ranking);
1185 free_evalue_refs(t);
1186 free(t);
1187 }
1188
1189 Q->next = next;
1190 }
1191
1192 Domain_Free(RD);
1193
1194 return ranking;
1195 }
1196
1197 Enumeration *barvinok_lexsmaller(Polyhedron *P, Polyhedron *D, unsigned dim,
1198 Polyhedron *C, unsigned MaxRays)
1199 {
1200 evalue *EP = barvinok_lexsmaller_ev(P, D, dim, C, MaxRays);
1201
1202 return partition2enumeration(EP);
1203 }
1204 #endif
1205
1206 /* "align" matrix to have nrows by inserting
1207 * the necessary number of rows and an equal number of columns in front
1208 */
1209 Matrix *align_matrix(Matrix *M, int nrows)
1210 {
1211 int i;
1212 int newrows = nrows - M->NbRows;
1213 Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
1214 for (i = 0; i < newrows; ++i)
1215 value_set_si(M2->p[i][i], 1);
1216 for (i = 0; i < M->NbRows; ++i)
1217 Vector_Copy(M->p[i], M2->p[newrows+i]+newrows, M->NbColumns);
1218 return M2;
1219 }
1220
1221 static void print_varlist(FILE *out, int n, char **names)
1222 {
1223 int i;
1224 fprintf(out, "[");
1225 for (i = 0; i < n; ++i) {
1226 if (i)
1227 fprintf(out, ",");
1228 fprintf(out, "%s", names[i]);
1229 }
1230 fprintf(out, "]");
1231 }
1232
1233 static void print_term(FILE *out, Value v, int pos, int dim, int nparam,
1234 char **iter_names, char **param_names, int *first)
1235 {
1236 if (value_zero_p(v)) {
1237 if (first && *first && pos >= dim + nparam)
1238 fprintf(out, "0");
1239 return;
1240 }
1241
1242 if (first) {
1243 if (!*first && value_pos_p(v))
1244 fprintf(out, "+");
1245 *first = 0;
1246 }
1247 if (pos < dim + nparam) {
1248 if (value_mone_p(v))
1249 fprintf(out, "-");
1250 else if (!value_one_p(v))
1251 value_print(out, VALUE_FMT, v);
1252 if (pos < dim)
1253 fprintf(out, "%s", iter_names[pos]);
1254 else
1255 fprintf(out, "%s", param_names[pos-dim]);
1256 } else
1257 value_print(out, VALUE_FMT, v);
1258 }
1259
1260 char **util_generate_names(int n, char *prefix)
1261 {
1262 int i;
1263 int len = (prefix ? strlen(prefix) : 0) + 10;
1264 char **names = ALLOCN(char*, n);
1265 if (!names) {
1266 fprintf(stderr, "ERROR: memory overflow.\n");
1267 exit(1);
1268 }
1269 for (i = 0; i < n; ++i) {
1270 names[i] = ALLOCN(char, len);
1271 if (!names[i]) {
1272 fprintf(stderr, "ERROR: memory overflow.\n");
1273 exit(1);
1274 }
1275 if (!prefix)
1276 snprintf(names[i], len, "%d", i);
1277 else
1278 snprintf(names[i], len, "%s%d", prefix, i);
1279 }
1280
1281 return names;
1282 }
1283
1284 void util_free_names(int n, char **names)
1285 {
1286 int i;
1287 for (i = 0; i < n; ++i)
1288 free(names[i]);
1289 free(names);
1290 }
1291
1292 void Polyhedron_pprint(FILE *out, Polyhedron *P, int dim, int nparam,
1293 char **iter_names, char **param_names)
1294 {
1295 int i, j;
1296 Value tmp;
1297
1298 assert(dim + nparam == P->Dimension);
1299
1300 value_init(tmp);
1301
1302 fprintf(out, "{ ");
1303 if (nparam) {
1304 print_varlist(out, nparam, param_names);
1305 fprintf(out, " -> ");
1306 }
1307 print_varlist(out, dim, iter_names);
1308 fprintf(out, " : ");
1309
1310 if (emptyQ2(P))
1311 fprintf(out, "FALSE");
1312 else for (i = 0; i < P->NbConstraints; ++i) {
1313 int first = 1;
1314 int v = First_Non_Zero(P->Constraint[i]+1, P->Dimension);
1315 if (v == -1 && value_pos_p(P->Constraint[i][0]))
1316 continue;
1317 if (i)
1318 fprintf(out, " && ");
1319 if (v == -1 && value_notzero_p(P->Constraint[i][1+P->Dimension]))
1320 fprintf(out, "FALSE");
1321 else if (value_pos_p(P->Constraint[i][v+1])) {
1322 print_term(out, P->Constraint[i][v+1], v, dim, nparam,
1323 iter_names, param_names, NULL);
1324 if (value_zero_p(P->Constraint[i][0]))
1325 fprintf(out, " = ");
1326 else
1327 fprintf(out, " >= ");
1328 for (j = v+1; j <= dim+nparam; ++j) {
1329 value_oppose(tmp, P->Constraint[i][1+j]);
1330 print_term(out, tmp, j, dim, nparam,
1331 iter_names, param_names, &first);
1332 }
1333 } else {
1334 value_oppose(tmp, P->Constraint[i][1+v]);
1335 print_term(out, tmp, v, dim, nparam,
1336 iter_names, param_names, NULL);
1337 fprintf(out, " <= ");
1338 for (j = v+1; j <= dim+nparam; ++j)
1339 print_term(out, P->Constraint[i][1+j], j, dim, nparam,
1340 iter_names, param_names, &first);
1341 }
1342 }
1343
1344 fprintf(out, " }\n");
1345
1346 value_clear(tmp);
1347 }
1348
1349 /* Construct a cone over P with P placed at x_d = 1, with
1350 * x_d the coordinate of an extra dimension
1351 *
1352 * It's probably a mistake to depend so much on the internal
1353 * representation. We should probably simply compute the
1354 * vertices/facets first.
1355 */
1356 Polyhedron *Cone_over_Polyhedron(Polyhedron *P)
1357 {
1358 unsigned NbConstraints = 0;
1359 unsigned NbRays = 0;
1360 Polyhedron *C;
1361 int i;
1362
1363 if (POL_HAS(P, POL_INEQUALITIES))
1364 NbConstraints = P->NbConstraints + 1;
1365 if (POL_HAS(P, POL_POINTS))
1366 NbRays = P->NbRays + 1;
1367
1368 C = Polyhedron_Alloc(P->Dimension+1, NbConstraints, NbRays);
1369 if (POL_HAS(P, POL_INEQUALITIES)) {
1370 C->NbEq = P->NbEq;
1371 for (i = 0; i < P->NbConstraints; ++i)
1372 Vector_Copy(P->Constraint[i], C->Constraint[i], P->Dimension+2);
1373 /* n >= 0 */
1374 value_set_si(C->Constraint[P->NbConstraints][0], 1);
1375 value_set_si(C->Constraint[P->NbConstraints][1+P->Dimension], 1);
1376 }
1377 if (POL_HAS(P, POL_POINTS)) {
1378 C->NbBid = P->NbBid;
1379 for (i = 0; i < P->NbRays; ++i)
1380 Vector_Copy(P->Ray[i], C->Ray[i], P->Dimension+2);
1381 /* vertex 0 */
1382 value_set_si(C->Ray[P->NbRays][0], 1);
1383 value_set_si(C->Ray[P->NbRays][1+C->Dimension], 1);
1384 }
1385 POL_SET(C, POL_VALID);
1386 if (POL_HAS(P, POL_INEQUALITIES))
1387 POL_SET(C, POL_INEQUALITIES);
1388 if (POL_HAS(P, POL_POINTS))
1389 POL_SET(C, POL_POINTS);
1390 if (POL_HAS(P, POL_VERTICES))
1391 POL_SET(C, POL_VERTICES);
1392 return C;
1393 }
1394
1395 const char *barvinok_version(void)
1396 {
1397 return
1398 "barvinok " VERSION " (" GIT_HEAD_ID ")\n"
1399 #ifdef USE_MODULO
1400 " +MODULO"
1401 #else
1402 " -MODULO"
1403 #endif
1404 #ifdef USE_INCREMENTAL_BF
1405 " INCREMENTAL=BF"
1406 #elif defined USE_INCREMENTAL_DF
1407 " INCREMENTAL=DF"
1408 #else
1409 " -INCREMENTAL"
1410 #endif
1411 "\n"
1412 #ifdef HAVE_CORRECT_VERTICES
1413 " +CORRECT_VERTICES"
1414 #else
1415 " -CORRECT_VERTICES"
1416 #endif
1417 #ifdef HAVE_PIPLIB
1418 " +PIPLIB"
1419 #else
1420 " -PIPLIB"
1421 #endif
1422 "\n"
1423 ;
1424 }
lattice point enumeration library