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