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