?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
regression test
[barvinok.git] / util.c
blob1d3f6d38914412ac4bb680e8b2c3011be53567eb
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 void manual_count(Polyhedron *P, Value* result)
11 {
12 Polyhedron *U = Universe_Polyhedron(0);
13 Enumeration *en = Polyhedron_Enumerate(P,U,1024,NULL);
14 Value *v = compute_poly(en,NULL);
15 value_assign(*result, *v);
16 value_clear(*v);
17 free(v);
18 Enumeration_Free(en);
19 Polyhedron_Free(U);
20 }
21
22 #include "ev_operations.h"
23 #include <util.h>
24 #include <barvinok.h>
25
26 /* Return random value between 0 and max-1 inclusive
27 */
28 int random_int(int max) {
29 return (int) (((double)(max))*rand()/(RAND_MAX+1.0));
30 }
31
32 /* Inplace polarization
33 */
34 void Polyhedron_Polarize(Polyhedron *P)
35 {
36 unsigned NbRows = P->NbConstraints + P->NbRays;
37 int i;
38 Value **q;
39
40 q = (Value **)malloc(NbRows * sizeof(Value *));
41 assert(q);
42 for (i = 0; i < P->NbRays; ++i)
43 q[i] = P->Ray[i];
44 for (; i < NbRows; ++i)
45 q[i] = P->Constraint[i-P->NbRays];
46 P->NbConstraints = NbRows - P->NbConstraints;
47 P->NbRays = NbRows - P->NbRays;
48 free(P->Constraint);
49 P->Constraint = q;
50 P->Ray = q + P->NbConstraints;
51 }
52
53 /*
54 * Rather general polar
55 * We can optimize it significantly if we assume that
56 * P includes zero
57 *
58 * Also, we calculate the polar as defined in Schrijver
59 * The opposite should probably work as well and would
60 * eliminate the need for multiplying by -1
61 */
62 Polyhedron* Polyhedron_Polar(Polyhedron *P, unsigned NbMaxRays)
63 {
64 int i;
65 Value mone;
66 unsigned dim = P->Dimension + 2;
67 Matrix *M = Matrix_Alloc(P->NbRays, dim);
68
69 assert(M);
70 value_init(mone);
71 value_set_si(mone, -1);
72 for (i = 0; i < P->NbRays; ++i) {
73 Vector_Scale(P->Ray[i], M->p[i], mone, dim);
74 value_multiply(M->p[i][0], M->p[i][0], mone);
75 value_multiply(M->p[i][dim-1], M->p[i][dim-1], mone);
76 }
77 P = Constraints2Polyhedron(M, NbMaxRays);
78 assert(P);
79 Matrix_Free(M);
80 value_clear(mone);
81 return P;
82 }
83
84 /*
85 * Returns the supporting cone of P at the vertex with index v
86 */
87 Polyhedron* supporting_cone(Polyhedron *P, int v)
88 {
89 Matrix *M;
90 Value tmp;
91 int i, n, j;
92 unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
93 unsigned dim = P->Dimension + 2;
94
95 assert(v >=0 && v < P->NbRays);
96 assert(value_pos_p(P->Ray[v][dim-1]));
97 assert(supporting);
98
99 value_init(tmp);
100 for (i = 0, n = 0; i < P->NbConstraints; ++i) {
101 Inner_Product(P->Constraint[i] + 1, P->Ray[v] + 1, dim - 1, &tmp);
102 if ((supporting[i] = value_zero_p(tmp)))
103 ++n;
104 }
105 assert(n >= dim - 2);
106 value_clear(tmp);
107 M = Matrix_Alloc(n, dim);
108 assert(M);
109 for (i = 0, j = 0; i < P->NbConstraints; ++i)
110 if (supporting[i]) {
111 value_set_si(M->p[j][dim-1], 0);
112 Vector_Copy(P->Constraint[i], M->p[j++], dim-1);
113 }
114 free(supporting);
115 P = Constraints2Polyhedron(M, P->NbRays+1);
116 assert(P);
117 Matrix_Free(M);
118 return P;
119 }
120
121 void value_lcm(Value i, Value j, Value* lcm)
122 {
123 Value aux;
124 value_init(aux);
125 value_multiply(aux,i,j);
126 Gcd(i,j,lcm);
127 value_division(*lcm,aux,*lcm);
128 value_clear(aux);
129 }
130
131 Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v)
132 {
133 Matrix *M;
134 Value lcm, tmp, tmp2;
135 unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
136 unsigned dim = P->Dimension + 2;
137 unsigned nparam = v->Vertex->NbColumns - 2;
138 unsigned nvar = dim - nparam - 2;
139 int i, n, j;
140 Vector *row;
141
142 assert(supporting);
143 row = Vector_Alloc(nparam+1);
144 assert(row);
145 value_init(lcm);
146 value_init(tmp);
147 value_init(tmp2);
148 value_set_si(lcm, 1);
149 for (i = 0, n = 0; i < P->NbConstraints; ++i) {
150 Vector_Set(row->p, 0, nparam+1);
151 for (j = 0 ; j < nvar; ++j) {
152 value_set_si(tmp, 1);
153 value_assign(tmp2, P->Constraint[i][j+1]);
154 if (value_ne(lcm, v->Vertex->p[j][nparam+1])) {
155 value_assign(tmp, lcm);
156 value_lcm(lcm, v->Vertex->p[j][nparam+1], &lcm);
157 value_division(tmp, lcm, tmp);
158 value_multiply(tmp2, tmp2, lcm);
159 value_division(tmp2, tmp2, v->Vertex->p[j][nparam+1]);
160 }
161 Vector_Combine(row->p, v->Vertex->p[j], row->p,
162 tmp, tmp2, nparam+1);
163 }
164 value_set_si(tmp, 1);
165 Vector_Combine(row->p, P->Constraint[i]+1+nvar, row->p, tmp, lcm, nparam+1);
166 for (j = 0; j < nparam+1; ++j)
167 if (value_notzero_p(row->p[j]))
168 break;
169 if ((supporting[i] = (j == nparam + 1)))
170 ++n;
171 }
172 assert(n >= nvar);
173 value_clear(tmp);
174 value_clear(tmp2);
175 value_clear(lcm);
176 Vector_Free(row);
177 M = Matrix_Alloc(n, nvar+2);
178 assert(M);
179 for (i = 0, j = 0; i < P->NbConstraints; ++i)
180 if (supporting[i]) {
181 value_set_si(M->p[j][nvar+1], 0);
182 Vector_Copy(P->Constraint[i], M->p[j++], nvar+1);
183 }
184 free(supporting);
185 P = Constraints2Polyhedron(M, P->NbRays+1);
186 assert(P);
187 Matrix_Free(M);
188 return P;
189 }
190
191 Polyhedron* triangularize_cone(Polyhedron *P, unsigned NbMaxCons)
192 {
193 const static int MAX_TRY=10;
194 int i, j, r, n, t;
195 Value tmp;
196 unsigned dim = P->Dimension;
197 Matrix *M = Matrix_Alloc(P->NbRays+1, dim+3);
198 Matrix *M2, *M3;
199 Polyhedron *L, *R, *T;
200 assert(P->NbEq == 0);
201
202 R = NULL;
203 value_init(tmp);
204
205 Vector_Set(M->p[0]+1, 0, dim+1);
206 value_set_si(M->p[0][0], 1);
207 value_set_si(M->p[0][dim+2], 1);
208 Vector_Set(M->p[P->NbRays]+1, 0, dim+2);
209 value_set_si(M->p[P->NbRays][0], 1);
210 value_set_si(M->p[P->NbRays][dim+1], 1);
211
212 for (i = 0, r = 1; i < P->NbRays; ++i) {
213 if (value_notzero_p(P->Ray[i][dim+1]))
214 continue;
215 Vector_Copy(P->Ray[i], M->p[r], dim+1);
216 Inner_Product(M->p[r]+1, M->p[r]+1, dim, &tmp);
217 value_assign(M->p[r][dim+1], tmp);
218 value_set_si(M->p[r][dim+2], 0);
219 ++r;
220 }
221
222 M3 = Matrix_Copy(M);
223 L = Rays2Polyhedron(M3, NbMaxCons);
224 Matrix_Free(M3);
225
226 M2 = Matrix_Alloc(dim+1, dim+2);
227
228 t = 1;
229 if (0) {
230 try_again:
231 /* Usually R should still be 0 */
232 Domain_Free(R);
233 Polyhedron_Free(L);
234 for (r = 1; r < P->NbRays; ++r) {
235 value_set_si(M->p[r][dim+1], random_int((t+1)*dim)+1);
236 }
237 M3 = Matrix_Copy(M);
238 L = Rays2Polyhedron(M3, NbMaxCons);
239 Matrix_Free(M3);
240 ++t;
241 }
242 assert(t <= MAX_TRY);
243
244 R = NULL;
245 n = 0;
246
247 for (i = 0; i < L->NbConstraints; ++i) {
248 if (value_negz_p(L->Constraint[i][dim+1]))
249 continue;
250 if (value_notzero_p(L->Constraint[i][dim+2]))
251 continue;
252 for (j = 1, r = 1; j < M->NbRows; ++j) {
253 Inner_Product(M->p[j]+1, L->Constraint[i]+1, dim+1, &tmp);
254 if (value_notzero_p(tmp))
255 continue;
256 if (r > dim)
257 goto try_again;
258 Vector_Copy(M->p[j]+1, M2->p[r]+1, dim);
259 value_set_si(M2->p[r][0], 1);
260 value_set_si(M2->p[r][dim+1], 0);
261 ++r;
262 }
263 assert(r == dim+1);
264 Vector_Set(M2->p[0]+1, 0, dim);
265 value_set_si(M2->p[0][0], 1);
266 value_set_si(M2->p[0][dim+1], 1);
267 T = Rays2Polyhedron(M2, P->NbConstraints+1);
268 T->next = R;
269 R = T;
270 ++n;
271 }
272 Matrix_Free(M2);
273
274 Polyhedron_Free(L);
275 value_clear(tmp);
276 Matrix_Free(M);
277
278 return R;
279 }
280
281 void check_triangulization(Polyhedron *P, Polyhedron *T)
282 {
283 Polyhedron *C, *D, *E, *F, *G, *U;
284 for (C = T; C; C = C->next) {
285 if (C == T)
286 U = C;
287 else
288 U = DomainConvex(DomainUnion(U, C, 100), 100);
289 for (D = C->next; D; D = D->next) {
290 F = C->next;
291 G = D->next;
292 C->next = NULL;
293 D->next = NULL;
294 E = DomainIntersection(C, D, 600);
295 assert(E->NbRays == 0 || E->NbEq >= 1);
296 Polyhedron_Free(E);
297 C->next = F;
298 D->next = G;
299 }
300 }
301 assert(PolyhedronIncludes(U, P));
302 assert(PolyhedronIncludes(P, U));
303 }
304
305 void Euclid(Value a, Value b, Value *x, Value *y, Value *g)
306 {
307 Value c, d, e, f, tmp;
308
309 value_init(c);
310 value_init(d);
311 value_init(e);
312 value_init(f);
313 value_init(tmp);
314 value_absolute(c, a);
315 value_absolute(d, b);
316 value_set_si(e, 1);
317 value_set_si(f, 0);
318 while(value_pos_p(d)) {
319 value_division(tmp, c, d);
320 value_multiply(tmp, tmp, f);
321 value_substract(e, e, tmp);
322 value_division(tmp, c, d);
323 value_multiply(tmp, tmp, d);
324 value_substract(c, c, tmp);
325 value_swap(c, d);
326 value_swap(e, f);
327 }
328 value_assign(*g, c);
329 if (value_zero_p(a))
330 value_set_si(*x, 0);
331 else if (value_pos_p(a))
332 value_assign(*x, e);
333 else value_oppose(*x, e);
334 if (value_zero_p(b))
335 value_set_si(*y, 0);
336 else {
337 value_multiply(tmp, a, *x);
338 value_substract(tmp, c, tmp);
339 value_division(*y, tmp, b);
340 }
341 value_clear(c);
342 value_clear(d);
343 value_clear(e);
344 value_clear(f);
345 value_clear(tmp);
346 }
347
348 Matrix * unimodular_complete(Vector *row)
349 {
350 Value g, b, c, old, tmp;
351 Matrix *m;
352 unsigned i, j;
353
354 value_init(b);
355 value_init(c);
356 value_init(g);
357 value_init(old);
358 value_init(tmp);
359 m = Matrix_Alloc(row->Size, row->Size);
360 for (j = 0; j < row->Size; ++j) {
361 value_assign(m->p[0][j], row->p[j]);
362 }
363 value_assign(g, row->p[0]);
364 for (i = 1; value_zero_p(g) && i < row->Size; ++i) {
365 for (j = 0; j < row->Size; ++j) {
366 if (j == i-1)
367 value_set_si(m->p[i][j], 1);
368 else
369 value_set_si(m->p[i][j], 0);
370 }
371 value_assign(g, row->p[i]);
372 }
373 for (; i < row->Size; ++i) {
374 value_assign(old, g);
375 Euclid(old, row->p[i], &c, &b, &g);
376 value_oppose(b, b);
377 for (j = 0; j < row->Size; ++j) {
378 if (j < i) {
379 value_multiply(tmp, row->p[j], b);
380 value_division(m->p[i][j], tmp, old);
381 } else if (j == i)
382 value_assign(m->p[i][j], c);
383 else
384 value_set_si(m->p[i][j], 0);
385 }
386 }
387 value_clear(b);
388 value_clear(c);
389 value_clear(g);
390 value_clear(old);
391 value_clear(tmp);
392 return m;
393 }
394
395 /*
396 * Returns a full-dimensional polyhedron with the same number
397 * of integer points as P
398 */
399 Polyhedron *remove_equalities(Polyhedron *P)
400 {
401 Value g;
402 Vector *v;
403 Polyhedron *p = Polyhedron_Copy(P), *q;
404 unsigned dim = p->Dimension;
405 Matrix *m1, *m2;
406 int i;
407
408 value_init(g);
409 while (p->NbEq > 0) {
410 assert(dim > 0);
411 Vector_Gcd(p->Constraint[0]+1, dim+1, &g);
412 Vector_AntiScale(p->Constraint[0]+1, p->Constraint[0]+1, g, dim+1);
413 Vector_Gcd(p->Constraint[0]+1, dim, &g);
414 if (value_notone_p(g) && value_notmone_p(g)) {
415 Polyhedron_Free(p);
416 p = Empty_Polyhedron(0);
417 break;
418 }
419 v = Vector_Alloc(dim);
420 Vector_Copy(p->Constraint[0]+1, v->p, dim);
421 m1 = unimodular_complete(v);
422 m2 = Matrix_Alloc(dim, dim+1);
423 for (i = 0; i < dim-1 ; ++i) {
424 Vector_Copy(m1->p[i+1], m2->p[i], dim);
425 value_set_si(m2->p[i][dim], 0);
426 }
427 Vector_Set(m2->p[dim-1], 0, dim);
428 value_set_si(m2->p[dim-1][dim], 1);
429 q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
430 Vector_Free(v);
431 Matrix_Free(m1);
432 Matrix_Free(m2);
433 Polyhedron_Free(p);
434 p = q;
435 --dim;
436 }
437 value_clear(g);
438 return p;
439 }
440
441 /*
442 * Returns a full-dimensional polyhedron with the same number
443 * of integer points as P
444 * nvar specifies the number of variables
445 * The remaining dimensions are assumed to be parameters
446 * Destroys P
447 * factor is NbEq x (nparam+2) matrix, containing stride constraints
448 * on the parameters; column nparam is the constant;
449 * column nparam+1 is the stride
450 */
451 Polyhedron *remove_equalities_p(Polyhedron *P, unsigned nvar, Matrix **factor)
452 {
453 Value g;
454 Vector *v;
455 Polyhedron *p = P, *q;
456 unsigned dim = p->Dimension;
457 Matrix *m1, *m2, *f;
458 int i, j, skip;
459
460 value_init(g);
461 f = Matrix_Alloc(p->NbEq, dim-nvar+2);
462 j = 0;
463 *factor = f;
464 skip = 0;
465 while (nvar > 0 && p->NbEq - skip > 0) {
466 assert(dim > 0);
467
468 while (value_zero_p(p->Constraint[skip][0]) &&
469 First_Non_Zero(p->Constraint[skip]+1, nvar) == -1)
470 ++skip;
471 if (p->NbEq == skip)
472 break;
473
474 Vector_Gcd(p->Constraint[skip]+1, dim+1, &g);
475 Vector_AntiScale(p->Constraint[skip]+1, p->Constraint[skip]+1, g, dim+1);
476 Vector_Gcd(p->Constraint[skip]+1, nvar, &g);
477 Vector_Copy(p->Constraint[skip]+1+nvar, f->p[j], dim-nvar+1);
478 value_assign(f->p[j][dim-nvar+1], g);
479 v = Vector_Alloc(dim);
480 Vector_AntiScale(p->Constraint[skip]+1, v->p, g, nvar);
481 Vector_Set(v->p+nvar, 0, dim-nvar);
482 m1 = unimodular_complete(v);
483 m2 = Matrix_Alloc(dim, dim+1);
484 for (i = 0; i < dim-1 ; ++i) {
485 Vector_Copy(m1->p[i+1], m2->p[i], dim);
486 value_set_si(m2->p[i][dim], 0);
487 }
488 Vector_Set(m2->p[dim-1], 0, dim);
489 value_set_si(m2->p[dim-1][dim], 1);
490 q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
491 Vector_Free(v);
492 Matrix_Free(m1);
493 Matrix_Free(m2);
494 Polyhedron_Free(p);
495 p = q;
496 --dim;
497 --nvar;
498 ++j;
499 }
500 value_clear(g);
501 return p;
502 }
503
504 struct single {
505 int nr;
506 int pos[2];
507 };
508
509 static void free_singles(int **singles, int dim)
510 {
511 int i;
512 for (i = 0; i < dim; ++i)
513 free(singles[i]);
514 free(singles);
515 }
516
517 static int **find_singles(Polyhedron *P, int dim, int max, int *nsingle)
518 {
519 int i, j, prev;
520 int **singles = (int **) malloc(dim * sizeof(int *));
521 assert(singles);
522
523 for (i = 0; i < dim; ++i) {
524 singles[i] = (int *) malloc((max + 1) *sizeof(int));
525 singles[i][0] = 0;
526 }
527
528 for (i = 0; i < P->NbConstraints; ++i) {
529 for (j = 0, prev = -1; j < dim; ++j) {
530 if (value_notzero_p(P->Constraint[i][j+1])) {
531 if (prev == -1)
532 prev = j;
533 else {
534 if (prev != -2)
535 singles[prev][0] = -1;
536 singles[j][0] = -1;
537 prev = -2;
538 }
539 }
540 }
541 if (prev >= 0 && singles[prev][0] >= 0)
542 singles[prev][++singles[prev][0]] = i;
543 }
544 *nsingle = 0;
545 for (j = 0; j < dim; ++j)
546 if (singles[j][0] > 0)
547 ++*nsingle;
548 if (!*nsingle) {
549 free_singles(singles, dim);
550 singles = 0;
551 }
552 return singles;
553 }
554
555 /*
556 * The number of points in P is equal to factor time
557 * the number of points in the polyhedron returned.
558 * The return value is zero if no reduction can be found.
559 */
560 Polyhedron* Polyhedron_Reduce(Polyhedron *P, Value* factor)
561 {
562 int i, j, nsingle, k, p;
563 unsigned dim = P->Dimension;
564 int **singles;
565
566 value_set_si(*factor, 1);
567
568 assert (P->NbEq == 0);
569
570 singles = find_singles(P, dim, 2, &nsingle);
571
572 if (nsingle == 0)
573 return 0;
574
575 {
576 Value tmp, pos, neg;
577 Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
578
579 value_init(tmp);
580 value_init(pos);
581 value_init(neg);
582
583 for (i = 0, j = 0; i < dim; ++i) {
584 if (singles[i][0] != 2)
585 value_set_si(m->p[j++][i], 1);
586 else {
587 for (k = 0; k <= 1; ++k) {
588 p = singles[i][1+k];
589 value_oppose(tmp, P->Constraint[p][dim+1]);
590 if (value_pos_p(P->Constraint[p][i+1]))
591 mpz_cdiv_q(pos, tmp, P->Constraint[p][i+1]);
592 else
593 mpz_fdiv_q(neg, tmp, P->Constraint[p][i+1]);
594 }
595 value_substract(tmp, neg, pos);
596 value_increment(tmp, tmp);
597 value_multiply(*factor, *factor, tmp);
598 }
599 }
600 value_set_si(m->p[dim-nsingle][dim], 1);
601 P = Polyhedron_Image(P, m, P->NbConstraints);
602 Matrix_Free(m);
603 free_singles(singles, dim);
604
605 value_clear(tmp);
606 value_clear(pos);
607 value_clear(neg);
608 }
609
610 return P;
611 }
612
613 /*
614 * Replaces constraint a x >= c by x >= ceil(c/a)
615 * where "a" is a common factor in the coefficients
616 * old is the constraint; v points to an initialized
617 * value that this procedure can use.
618 * Return non-zero if something changed.
619 * Result is places in new.
620 */
621 int ConstraintSimplify(Value *old, Value *new, int len, Value* v)
622 {
623 Vector_Gcd(old+1, len - 2, v);
624
625 if (value_one_p(*v))
626 return 0;
627
628 Vector_AntiScale(old+1, new+1, *v, len-2);
629 mpz_fdiv_q(new[len-1], old[len-1], *v);
630 return 1;
631 }
632
633 /*
634 * Project on final dim dimensions
635 */
636 static Polyhedron* Polyhedron_Project(Polyhedron *P, int dim)
637 {
638 int i;
639 int remove = P->Dimension - dim;
640
641 if (P->Dimension == dim)
642 return Polyhedron_Copy(P);
643
644 Matrix *T = Matrix_Alloc(dim+1, P->Dimension+1);
645 for (i = 0; i < dim+1; ++i)
646 value_set_si(T->p[i][i+remove], 1);
647 Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
648 Matrix_Free(T);
649 return I;
650 }
651
652 struct section { Polyhedron * D; evalue E; };
653
654 static Polyhedron* ParamPolyhedron_Reduce_mod(Polyhedron *P, unsigned nvar,
655 evalue* factor)
656 {
657 int nsingle;
658 int **singles;
659 unsigned dim = P->Dimension;
660
661 singles = find_singles(P, nvar, P->NbConstraints, &nsingle);
662
663 if (nsingle == 0)
664 return 0;
665
666 {
667 Polyhedron *C, *T;
668 int i, j, p, n;
669 Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
670 Value tmp, g;
671 evalue mone;
672 value_init(mone.d);
673 evalue_set_si(&mone, -1, 1);
674 C = Polyhedron_Project(P, dim-nvar);
675
676 value_init(tmp);
677 value_init(g);
678
679 for (i = 0, j = 0; i < dim; ++i) {
680 if (i >= nvar || singles[i][0] < 2)
681 value_set_si(m->p[j++][i], 1);
682 else {
683 struct section *s;
684 Matrix *M, *M2;
685 int nd = 0;
686 int k, l, k2, l2, q;
687 evalue *L, *U;
688 evalue F;
689 /* put those with positive coefficients first; number: p */
690 for (p = 0, n = singles[i][0]-1; p <= n; ) {
691 while (value_pos_p(P->Constraint[singles[i][1+p]][i+1]))
692 ++p;
693 while (value_neg_p(P->Constraint[singles[i][1+n]][i+1]))
694 --n;
695 if (p < n) {
696 int t = singles[i][1+p];
697 singles[i][1+p] = singles[i][1+n];
698 singles[i][1+n] = t;
699 ++p;
700 --n;
701 }
702 }
703 n = singles[i][0]-p;
704 assert (p >= 1 && n >= 1);
705 s = (struct section *) malloc(p * n * sizeof(struct section));
706 M = Matrix_Alloc((p-1) + (n-1), dim-nvar+2);
707 for (k = 0; k < p; ++k) {
708 for (k2 = 0; k2 < p; ++k2) {
709 if (k2 == k)
710 continue;
711 q = k2 - (k2 > k);
712 value_oppose(tmp, P->Constraint[singles[i][1+k2]][i+1]);
713 value_set_si(M->p[q][0], 1);
714 Vector_Combine(P->Constraint[singles[i][1+k]]+1+nvar,
715 P->Constraint[singles[i][1+k2]]+1+nvar,
716 M->p[q]+1,
717 tmp,
718 P->Constraint[singles[i][1+k]][i+1],
719 dim-nvar+1);
720 if (k2 > k)
721 value_decrement(M->p[q][dim-nvar+1],
722 M->p[q][dim-nvar+1]);
723 ConstraintSimplify(M->p[q], M->p[q],
724 dim-nvar+2, &g);
725 }
726 for (l = p; l < p+n; ++l) {
727 value_oppose(tmp, P->Constraint[singles[i][1+l]][i+1]);
728 for (l2 = p; l2 < p+n; ++l2) {
729 if (l2 == l)
730 continue;
731 q = l2-1 - (l2 > l);
732 value_set_si(M->p[q][0], 1);
733 Vector_Combine(P->Constraint[singles[i][1+l2]]+1+nvar,
734 P->Constraint[singles[i][1+l]]+1+nvar,
735 M->p[q]+1,
736 tmp,
737 P->Constraint[singles[i][1+l2]][i+1],
738 dim-nvar+1);
739 if (l2 > l)
740 value_decrement(M->p[q][dim-nvar+1],
741 M->p[q][dim-nvar+1]);
742 ConstraintSimplify(M->p[q], M->p[q],
743 dim-nvar+2, &g);
744 }
745 M2 = Matrix_Copy(M);
746 s[nd].D = Constraints2Polyhedron(M2, P->NbRays);
747 Matrix_Free(M2);
748 if (emptyQ(s[nd].D)) {
749 Polyhedron_Free(s[nd].D);
750 continue;
751 }
752 T = DomainIntersection(s[nd].D, C, 0);
753 L = bv_ceil3(P->Constraint[singles[i][1+k]]+1+nvar,
754 dim-nvar+1,
755 P->Constraint[singles[i][1+k]][i+1], T);
756 U = bv_ceil3(P->Constraint[singles[i][1+l]]+1+nvar,
757 dim-nvar+1,
758 P->Constraint[singles[i][1+l]][i+1], T);
759 Domain_Free(T);
760 eadd(L, U);
761 eadd(&mone, U);
762 emul(&mone, U);
763 s[nd].E = *U;
764 free_evalue_refs(L);
765 free(L);
766 free(U);
767 ++nd;
768 }
769 }
770
771 Matrix_Free(M);
772
773 value_init(F.d);
774 value_set_si(F.d, 0);
775 F.x.p = new_enode(partition, 2*nd, -1);
776 for (k = 0; k < nd; ++k) {
777 EVALUE_SET_DOMAIN(F.x.p->arr[2*k], s[k].D);
778 value_clear(F.x.p->arr[2*k+1].d);
779 F.x.p->arr[2*k+1] = s[k].E;
780 }
781 free(s);
782
783 emul(&F, factor);
784 free_evalue_refs(&F);
785 }
786 }
787 value_set_si(m->p[dim-nsingle][dim], 1);
788 P = Polyhedron_Image(P, m, P->NbConstraints);
789 Matrix_Free(m);
790 free_singles(singles, nvar);
791
792 value_clear(g);
793 value_clear(tmp);
794
795 free_evalue_refs(&mone);
796 Polyhedron_Free(C);
797 }
798
799 reduce_evalue(factor);
800
801 return P;
802 }
803
804 #ifdef USE_MODULO
805 Polyhedron* ParamPolyhedron_Reduce(Polyhedron *P, unsigned nvar,
806 evalue* factor)
807 {
808 return ParamPolyhedron_Reduce_mod(P, nvar, factor);
809 }
810 #else
811 Polyhedron* ParamPolyhedron_Reduce(Polyhedron *P, unsigned nvar,
812 evalue* factor)
813 {
814 Polyhedron *R;
815 evalue tmp;
816 value_init(tmp.d);
817 evalue_set_si(&tmp, 1, 1);
818 R = ParamPolyhedron_Reduce_mod(P, nvar, &tmp);
819 evalue_mod2table(&tmp, P->Dimension - nvar);
820 reduce_evalue(&tmp);
821 emul(&tmp, factor);
822 free_evalue_refs(&tmp);
823 return R;
824 }
825 #endif
826
827 Bool isIdentity(Matrix *M)
828 {
829 unsigned i, j;
830 if (M->NbRows != M->NbColumns)
831 return False;
832
833 for (i = 0;i < M->NbRows; i ++)
834 for (j = 0; j < M->NbColumns; j ++)
835 if (i == j) {
836 if(value_notone_p(M->p[i][j]))
837 return False;
838 } else {
839 if(value_notzero_p(M->p[i][j]))
840 return False;
841 }
842 return True;
843 }
844
845 void Param_Polyhedron_Print(FILE* DST, Param_Polyhedron *PP, char **param_names)
846 {
847 Param_Domain *P;
848 Param_Vertices *V;
849
850 for(P=PP->D;P;P=P->next) {
851
852 /* prints current val. dom. */
853 printf( "---------------------------------------\n" );
854 printf( "Domain :\n");
855 Print_Domain( stdout, P->Domain, param_names );
856
857 /* scan the vertices */
858 printf( "Vertices :\n");
859 FORALL_PVertex_in_ParamPolyhedron(V,P,PP) {
860
861 /* prints each vertex */
862 Print_Vertex( stdout, V->Vertex, param_names );
863 printf( "\n" );
864 }
865 END_FORALL_PVertex_in_ParamPolyhedron;
866 }
867 }
868
869 void Enumeration_Print(FILE *Dst, Enumeration *en, char **params)
870 {
871 for (; en; en = en->next) {
872 Print_Domain(Dst, en->ValidityDomain, params);
873 print_evalue(Dst, &en->EP, params);
874 }
875 }
876
877 void Enumeration_mod2table(Enumeration *en, unsigned nparam)
878 {
879 for (; en; en = en->next) {
880 evalue_mod2table(&en->EP, nparam);
881 reduce_evalue(&en->EP);
882 }
883 }
884
885 size_t Enumeration_size(Enumeration *en)
886 {
887 size_t s = 0;
888
889 for (; en; en = en->next) {
890 s += domain_size(en->ValidityDomain);
891 s += evalue_size(&en->EP);
892 }
893 return s;
894 }
895
896 void Free_ParamNames(char **params, int m)
897 {
898 while (--m >= 0)
899 free(params[m]);
900 free(params);
901 }
902
903 int DomainIncludes(Polyhedron *Pol1, Polyhedron *Pol2)
904 {
905 Polyhedron *P2;
906 for ( ; Pol1; Pol1 = Pol1->next) {
907 for (P2 = Pol2; P2; P2 = P2->next)
908 if (!PolyhedronIncludes(Pol1, P2))
909 break;
910 if (!P2)
911 return 1;
912 }
913 return 0;
914 }
915
916 static Polyhedron *p_simplify_constraints(Polyhedron *P, Vector *row,
917 Value *g, unsigned MaxRays)
918 {
919 Polyhedron *T, *R = P;
920 int len = P->Dimension+2;
921 int r;
922
923 /* Also look at equalities.
924 * If an equality can be "simplified" then there
925 * are no integer solutions anyway and the following loop
926 * will add a conflicting constraint
927 */
928 for (r = 0; r < R->NbConstraints; ++r) {
929 if (ConstraintSimplify(R->Constraint[r], row->p, len, g)) {
930 T = R;
931 R = AddConstraints(row->p, 1, R, MaxRays);
932 if (T != P)
933 Polyhedron_Free(T);
934 r = -1;
935 }
936 }
937 if (R != P)
938 Polyhedron_Free(P);
939 return R;
940 }
941
942 /*
943 * Replaces constraint a x >= c by x >= ceil(c/a)
944 * where "a" is a common factor in the coefficients
945 * Destroys P and returns a newly allocated Polyhedron
946 * or just returns P in case no changes were made
947 */
948 Polyhedron *DomainConstraintSimplify(Polyhedron *P, unsigned MaxRays)
949 {
950 Polyhedron **prev;
951 int len = P->Dimension+2;
952 Vector *row = Vector_Alloc(len);
953 value_set_si(row->p[0], 1);
954 Value g;
955 value_init(g);
956 Polyhedron *R = P, *N;
957
958 for (prev = &R; P; P = N) {
959 Polyhedron *T;
960 N = P->next;
961 T = p_simplify_constraints(P, row, &g, MaxRays);
962
963 if (emptyQ(T) && prev != &R) {
964 Polyhedron_Free(T);
965 *prev = NULL;
966 continue;
967 }
968
969 if (T != P)
970 T->next = N;
971 *prev = T;
972 prev = &T->next;
973 }
974
975 if (R->next && emptyQ(R)) {
976 N = R->next;
977 Polyhedron_Free(R);
978 R = N;
979 }
980
981 value_clear(g);
982 Vector_Free(row);
983 return R;
984 }
985
986 void line_minmax(Polyhedron *I, Value *min, Value *max)
987 {
988 int i;
989
990 if (I->NbEq >= 1) {
991 value_oppose(I->Constraint[0][2], I->Constraint[0][2]);
992 /* There should never be a remainder here */
993 if (value_pos_p(I->Constraint[0][1]))
994 mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
995 else
996 mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
997 value_assign(*max, *min);
998 } else for (i = 0; i < I->NbConstraints; ++i) {
999 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
1000 if (value_pos_p(I->Constraint[i][1]))
1001 mpz_cdiv_q(*min, I->Constraint[i][2], I->Constraint[i][1]);
1002 else
1003 mpz_fdiv_q(*max, I->Constraint[i][2], I->Constraint[i][1]);
1004 }
1005 Polyhedron_Free(I);
1006 }
lattice point enumeration library