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