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