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