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