?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
iscc: support transitive closure
[barvinok.git] / barvinok_e.cc
blob9292ad1e8b9c4c7223f559e5b3d0afb5148b5c9d
1 #include <assert.h>
2 #include <isl_set_polylib.h>
3 #include <barvinok/barvinok.h>
4 #include <barvinok/evalue.h>
5 #include <barvinok/util.h>
6 #include "param_util.h"
7 #include "piputil.h"
8 #include "reduce_domain.h"
9 #include "config.h"
10
11 #define ALLOC(type) (type*)malloc(sizeof(type))
12
13 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
14 {
15 int len = P->Dimension+2;
16 Polyhedron *T, *R = P;
17 Value g;
18 value_init(g);
19 Vector *row = Vector_Alloc(len);
20 value_set_si(row->p[0], 1);
21
22 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
23
24 Matrix *M = Matrix_Alloc(2, len-1);
25 value_set_si(M->p[1][len-2], 1);
26 for (int v = 0; v < P->Dimension; ++v) {
27 value_set_si(M->p[0][v], 1);
28 Polyhedron *I = Polyhedron_Image(R, M, 2+1);
29 value_set_si(M->p[0][v], 0);
30 for (int r = 0; r < I->NbConstraints; ++r) {
31 if (value_zero_p(I->Constraint[r][0]))
32 continue;
33 if (value_zero_p(I->Constraint[r][1]))
34 continue;
35 if (value_one_p(I->Constraint[r][1]))
36 continue;
37 if (value_mone_p(I->Constraint[r][1]))
38 continue;
39 value_absolute(g, I->Constraint[r][1]);
40 Vector_Set(row->p+1, 0, len-2);
41 value_division(row->p[1+v], I->Constraint[r][1], g);
42 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
43 T = R;
44 R = AddConstraints(row->p, 1, R, MaxRays);
45 if (T != P)
46 Polyhedron_Free(T);
47 }
48 Polyhedron_Free(I);
49 }
50 Matrix_Free(M);
51 Vector_Free(row);
52 value_clear(g);
53 return R;
54 }
55
56 /* Construct a constraint c from constraints l and u such that if
57 * if constraint c holds then for each value of the other variables
58 * there is at most one value of variable pos (position pos+1 in the constraints).
59 *
60 * Given a lower and an upper bound
61 * n_l v_i + <c_l,x> + c_l >= 0
62 * -n_u v_i + <c_u,x> + c_u >= 0
63 * the constructed constraint is
64 *
65 * -(n_l<c_u,x> + n_u<c_l,x>) + (-n_l c_u - n_u c_l + n_l n_u - 1)
66 *
67 * which is then simplified to remove the content of the non-constant coefficients
68 *
69 * len is the total length of the constraints.
70 * v is a temporary variable that can be used by this procedure
71 */
72 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
73 int len, Value *v)
74 {
75 value_oppose(*v, u[pos+1]);
76 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
77 value_multiply(*v, *v, l[pos+1]);
78 value_subtract(c[len-1], c[len-1], *v);
79 value_set_si(*v, -1);
80 Vector_Scale(c+1, c+1, *v, len-1);
81 value_decrement(c[len-1], c[len-1]);
82 ConstraintSimplify(c, c, len, v);
83 }
84
85 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
86 int len)
87 {
88 bool parallel;
89 Value g1;
90 Value g2;
91 value_init(g1);
92 value_init(g2);
93
94 Vector_Gcd(&l[1+pos], len, &g1);
95 Vector_Gcd(&u[1+pos], len, &g2);
96 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
97 parallel = First_Non_Zero(c+1, len) == -1;
98
99 value_clear(g1);
100 value_clear(g2);
101
102 return parallel;
103 }
104
105 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
106 int exist, int len, Value *v)
107 {
108 Value g;
109 value_init(g);
110
111 Vector_Gcd(&u[1+pos], exist, v);
112 Vector_Gcd(&l[1+pos], exist, &g);
113 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
114 value_multiply(*v, *v, g);
115 value_subtract(c[len-1], c[len-1], *v);
116 value_set_si(*v, -1);
117 Vector_Scale(c+1, c+1, *v, len-1);
118 value_decrement(c[len-1], c[len-1]);
119 ConstraintSimplify(c, c, len, v);
120
121 value_clear(g);
122 }
123
124 /* Turns a x + b >= 0 into a x + b <= -1
125 *
126 * len is the total length of the constraint.
127 * v is a temporary variable that can be used by this procedure
128 */
129 static void oppose_constraint(Value *c, int len, Value *v)
130 {
131 value_set_si(*v, -1);
132 Vector_Scale(c+1, c+1, *v, len-1);
133 value_decrement(c[len-1], c[len-1]);
134 }
135
136 /* Split polyhedron P into two polyhedra *pos and *neg, where
137 * existential variable i has at most one solution for each
138 * value of the other variables in *neg.
139 *
140 * The splitting is performed using constraints l and u.
141 *
142 * nvar: number of set variables
143 * row: temporary vector that can be used by this procedure
144 * f: temporary value that can be used by this procedure
145 */
146 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
147 int nvar, int MaxRays, Vector *row, Value& f,
148 Polyhedron **pos, Polyhedron **neg)
149 {
150 negative_test_constraint(P->Constraint[l], P->Constraint[u],
151 row->p, nvar+i, P->Dimension+2, &f);
152 *neg = AddConstraints(row->p, 1, P, MaxRays);
153 POL_ENSURE_VERTICES(*neg);
154
155 /* We found an independent, but useless constraint
156 * Maybe we should detect this earlier and not
157 * mark the variable as INDEPENDENT
158 */
159 if (emptyQ((*neg))) {
160 Polyhedron_Free(*neg);
161 return false;
162 }
163
164 oppose_constraint(row->p, P->Dimension+2, &f);
165 *pos = AddConstraints(row->p, 1, P, MaxRays);
166 POL_ENSURE_VERTICES(*pos);
167
168 if (emptyQ((*pos))) {
169 Polyhedron_Free(*neg);
170 Polyhedron_Free(*pos);
171 return false;
172 }
173
174 return true;
175 }
176
177 /*
178 * unimodularly transform P such that constraint r is transformed
179 * into a constraint that involves only a single (the first)
180 * existential variable
181 *
182 */
183 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
184 unsigned MaxRays)
185 {
186 Value g;
187 value_init(g);
188
189 Matrix *M = Matrix_Alloc(exist, exist);
190 Vector_Copy(P->Constraint[r]+1+nvar, M->p[0], exist);
191 Vector_Gcd(M->p[0], exist, &g);
192 if (value_notone_p(g))
193 Vector_AntiScale(M->p[0], M->p[0], g, exist);
194 value_clear(g);
195
196 int ok = unimodular_complete(M, 1);
197 assert(ok);
198 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
199 for (r = 0; r < nvar; ++r)
200 value_set_si(M2->p[r][r], 1);
201 for ( ; r < nvar+exist; ++r)
202 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
203 for ( ; r < P->Dimension+1; ++r)
204 value_set_si(M2->p[r][r], 1);
205 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
206
207 Matrix_Free(M2);
208 Matrix_Free(M);
209
210 return T;
211 }
212
213 /* Split polyhedron P into two polyhedra *pos and *neg, where
214 * existential variable i has at most one solution for each
215 * value of the other variables in *neg.
216 *
217 * If independent is set, then the two constraints on which the
218 * split will be performed need to be independent of the other
219 * existential variables.
220 *
221 * Return true if an appropriate split could be performed.
222 *
223 * nvar: number of set variables
224 * exist: number of existential variables
225 * row: temporary vector that can be used by this procedure
226 * f: temporary value that can be used by this procedure
227 */
228 static bool SplitOnVar(Polyhedron *P, int i,
229 int nvar, int exist, int MaxRays,
230 Vector *row, Value& f, bool independent,
231 Polyhedron **pos, Polyhedron **neg)
232 {
233 int j;
234
235 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
236 if (value_negz_p(P->Constraint[l][nvar+i+1]))
237 continue;
238
239 if (independent) {
240 for (j = 0; j < exist; ++j)
241 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
242 break;
243 if (j < exist)
244 continue;
245 }
246
247 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
248 if (value_posz_p(P->Constraint[u][nvar+i+1]))
249 continue;
250
251 if (independent) {
252 for (j = 0; j < exist; ++j)
253 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
254 break;
255 if (j < exist)
256 continue;
257 }
258
259 if (SplitOnConstraint(P, i, l, u, nvar, MaxRays, row, f, pos, neg)) {
260 if (independent) {
261 if (i != 0)
262 Polyhedron_ExchangeColumns(*neg, nvar+1, nvar+1+i);
263 }
264 return true;
265 }
266 }
267 }
268
269 return false;
270 }
271
272 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
273 int i, int l1, int l2,
274 Polyhedron **pos, Polyhedron **neg)
275 {
276 Value f;
277 value_init(f);
278 Vector *row = Vector_Alloc(P->Dimension+2);
279 value_set_si(row->p[0], 1);
280 value_oppose(f, P->Constraint[l1][nvar+i+1]);
281 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
282 row->p+1,
283 P->Constraint[l2][nvar+i+1], f,
284 P->Dimension+1);
285 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
286 *pos = AddConstraints(row->p, 1, P, 0);
287 POL_ENSURE_VERTICES(*pos);
288 value_set_si(f, -1);
289 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
290 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
291 *neg = AddConstraints(row->p, 1, P, 0);
292 POL_ENSURE_VERTICES(*neg);
293 Vector_Free(row);
294 value_clear(f);
295
296 return !emptyQ((*pos)) && !emptyQ((*neg));
297 }
298
299 static bool double_bound(Polyhedron *P, int nvar, int exist,
300 Polyhedron **pos, Polyhedron **neg)
301 {
302 for (int i = 0; i < exist; ++i) {
303 int l1, l2;
304 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
305 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
306 continue;
307 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
308 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
309 continue;
310 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
311 return true;
312 }
313 }
314 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
315 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
316 continue;
317 if (l1 < P->NbConstraints)
318 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
319 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
320 continue;
321 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
322 return true;
323 }
324 }
325 return false;
326 }
327 return false;
328 }
329
330 enum constraint {
331 ALL_POS = 1 << 0,
332 ONE_NEG = 1 << 1,
333 INDEPENDENT = 1 << 2,
334 ROT_NEG = 1 << 3
335 };
336
337 static evalue* enumerate_or(Polyhedron *D,
338 unsigned exist, unsigned nparam, barvinok_options *options)
339 {
340 #ifdef DEBUG_ER
341 fprintf(stderr, "\nER: Or\n");
342 #endif /* DEBUG_ER */
343
344 Polyhedron *N = D->next;
345 D->next = 0;
346 evalue *EP =
347 barvinok_enumerate_e_with_options(D, exist, nparam, options);
348 Polyhedron_Free(D);
349
350 for (D = N; D; D = N) {
351 N = D->next;
352 D->next = 0;
353
354 evalue *EN =
355 barvinok_enumerate_e_with_options(D, exist, nparam, options);
356
357 eor(EN, EP);
358 evalue_free(EN);
359 Polyhedron_Free(D);
360 }
361
362 reduce_evalue(EP);
363
364 return EP;
365 }
366
367 static evalue* enumerate_sum(Polyhedron *P,
368 unsigned exist, unsigned nparam, barvinok_options *options)
369 {
370 int nvar = P->Dimension - exist - nparam;
371 int toswap = nvar < exist ? nvar : exist;
372 for (int i = 0; i < toswap; ++i)
373 Polyhedron_ExchangeColumns(P, 1 + i, nvar+exist - i);
374 nparam += nvar;
375
376 #ifdef DEBUG_ER
377 fprintf(stderr, "\nER: Sum\n");
378 #endif /* DEBUG_ER */
379
380 evalue *EP = barvinok_enumerate_e_with_options(P, exist, nparam, options);
381
382 evalue_split_domains_into_orthants(EP, options->MaxRays);
383 reduce_evalue(EP);
384 evalue_range_reduction(EP);
385
386 evalue_frac2floor(EP);
387
388 evalue *sum = barvinok_summate(EP, nvar, options);
389
390 evalue_free(EP);
391 EP = sum;
392
393 evalue_range_reduction(EP);
394
395 return EP;
396 }
397
398 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
399 unsigned exist, unsigned nparam, barvinok_options *options)
400 {
401 int nvar = P->Dimension - exist - nparam;
402
403 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
404 for (int i = 0; i < exist; ++i)
405 value_set_si(M->p[i][nvar+i+1], 1);
406 Polyhedron *O = S;
407 S = DomainAddRays(S, M, options->MaxRays);
408 Polyhedron_Free(O);
409 Polyhedron *F = DomainAddRays(P, M, options->MaxRays);
410 Polyhedron *D = DomainDifference(F, S, options->MaxRays);
411 O = D;
412 D = Disjoint_Domain(D, 0, options->MaxRays);
413 Polyhedron_Free(F);
414 Domain_Free(O);
415 Matrix_Free(M);
416
417 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
418 for (int j = 0; j < nvar; ++j)
419 value_set_si(M->p[j][j], 1);
420 for (int j = 0; j < nparam+1; ++j)
421 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
422 Polyhedron *T = Polyhedron_Image(S, M, options->MaxRays);
423 evalue *EP = barvinok_enumerate_e_with_options(T, 0, nparam, options);
424 Polyhedron_Free(S);
425 Polyhedron_Free(T);
426 Matrix_Free(M);
427
428 for (Polyhedron *Q = D; Q; Q = Q->next) {
429 Polyhedron *N = Q->next;
430 Q->next = 0;
431 T = DomainIntersection(P, Q, options->MaxRays);
432 evalue *E = barvinok_enumerate_e_with_options(T, exist, nparam, options);
433 eadd(E, EP);
434 evalue_free(E);
435 Polyhedron_Free(T);
436 Q->next = N;
437 }
438 Domain_Free(D);
439 return EP;
440 }
441
442 static evalue* enumerate_sure(Polyhedron *P,
443 unsigned exist, unsigned nparam, barvinok_options *options)
444 {
445 int i;
446 Polyhedron *S = P;
447 int nvar = P->Dimension - exist - nparam;
448 Value lcm;
449 Value f;
450 value_init(lcm);
451 value_init(f);
452
453 for (i = 0; i < exist; ++i) {
454 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
455 int c = 0;
456 value_set_si(lcm, 1);
457 for (int j = 0; j < S->NbConstraints; ++j) {
458 if (value_negz_p(S->Constraint[j][1+nvar+i]))
459 continue;
460 if (value_one_p(S->Constraint[j][1+nvar+i]))
461 continue;
462 value_lcm(lcm, lcm, S->Constraint[j][1+nvar+i]);
463 }
464
465 for (int j = 0; j < S->NbConstraints; ++j) {
466 if (value_negz_p(S->Constraint[j][1+nvar+i]))
467 continue;
468 if (value_one_p(S->Constraint[j][1+nvar+i]))
469 continue;
470 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
471 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
472 value_subtract(M->p[c][S->Dimension+1],
473 M->p[c][S->Dimension+1],
474 lcm);
475 value_increment(M->p[c][S->Dimension+1],
476 M->p[c][S->Dimension+1]);
477 ++c;
478 }
479 Polyhedron *O = S;
480 S = AddConstraints(M->p[0], c, S, options->MaxRays);
481 if (O != P)
482 Polyhedron_Free(O);
483 Matrix_Free(M);
484 if (emptyQ(S)) {
485 Polyhedron_Free(S);
486 value_clear(lcm);
487 value_clear(f);
488 return 0;
489 }
490 }
491 value_clear(lcm);
492 value_clear(f);
493
494 #ifdef DEBUG_ER
495 fprintf(stderr, "\nER: Sure\n");
496 #endif /* DEBUG_ER */
497
498 return split_sure(P, S, exist, nparam, options);
499 }
500
501 static evalue* enumerate_sure2(Polyhedron *P,
502 unsigned exist, unsigned nparam, barvinok_options *options)
503 {
504 int nvar = P->Dimension - exist - nparam;
505 int r;
506 for (r = 0; r < P->NbRays; ++r)
507 if (value_one_p(P->Ray[r][0]) &&
508 value_one_p(P->Ray[r][P->Dimension+1]))
509 break;
510
511 if (r >= P->NbRays)
512 return 0;
513
514 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
515 for (int i = 0; i < nvar; ++i)
516 value_set_si(M->p[i][1+i], 1);
517 for (int i = 0; i < nparam; ++i)
518 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
519 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
520 value_set_si(M->p[nvar+nparam][0], 1);
521 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
522 Polyhedron * F = Rays2Polyhedron(M, options->MaxRays);
523 Matrix_Free(M);
524
525 Polyhedron *I = DomainIntersection(F, P, options->MaxRays);
526 Polyhedron_Free(F);
527
528 #ifdef DEBUG_ER
529 fprintf(stderr, "\nER: Sure2\n");
530 #endif /* DEBUG_ER */
531
532 return split_sure(P, I, exist, nparam, options);
533 }
534
535 static evalue* enumerate_cyclic(Polyhedron *P,
536 unsigned exist, unsigned nparam,
537 evalue * EP, int r, int p, unsigned MaxRays)
538 {
539 int nvar = P->Dimension - exist - nparam;
540
541 /* If EP in its fractional maps only contains references
542 * to the remainder parameter with appropriate coefficients
543 * then we could in principle avoid adding existentially
544 * quantified variables to the validity domains.
545 * We'd have to replace the remainder by m { p/m }
546 * and multiply with an appropriate factor that is one
547 * only in the appropriate range.
548 * This last multiplication can be avoided if EP
549 * has a single validity domain with no (further)
550 * constraints on the remainder parameter
551 */
552
553 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
554 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
555 for (int j = 0; j < nparam; ++j)
556 if (j != p)
557 value_set_si(CT->p[j][j], 1);
558 value_set_si(CT->p[p][nparam+1], 1);
559 value_set_si(CT->p[nparam][nparam+2], 1);
560 value_set_si(M->p[0][1+p], -1);
561 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
562 value_set_si(M->p[0][1+nparam+1], 1);
563 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
564 Matrix_Free(M);
565 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
566 Polyhedron_Free(CEq);
567 Matrix_Free(CT);
568
569 return EP;
570 }
571
572 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
573 {
574 if (value_notzero_p(EP->d))
575 return;
576
577 assert(EP->x.p->type == partition);
578 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
579 for (int i = 0; i < EP->x.p->size/2; ++i) {
580 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
581 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
582 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
583 Domain_Free(D);
584 }
585 }
586
587 static evalue* enumerate_line(Polyhedron *P,
588 unsigned exist, unsigned nparam, barvinok_options *options)
589 {
590 if (P->NbBid == 0)
591 return 0;
592
593 #ifdef DEBUG_ER
594 fprintf(stderr, "\nER: Line\n");
595 #endif /* DEBUG_ER */
596
597 int nvar = P->Dimension - exist - nparam;
598 int i, j;
599 for (i = 0; i < nparam; ++i)
600 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
601 break;
602 assert(i < nparam);
603 for (j = i+1; j < nparam; ++j)
604 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
605 break;
606 assert(j >= nparam); // for now
607
608 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
609 value_set_si(M->p[0][0], 1);
610 value_set_si(M->p[0][1+nvar+exist+i], 1);
611 value_set_si(M->p[1][0], 1);
612 value_set_si(M->p[1][1+nvar+exist+i], -1);
613 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
614 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
615 Polyhedron *S = AddConstraints(M->p[0], 2, P, options->MaxRays);
616 evalue *EP = barvinok_enumerate_e_with_options(S, exist, nparam, options);
617 Polyhedron_Free(S);
618 Matrix_Free(M);
619
620 return enumerate_cyclic(P, exist, nparam, EP, 0, i, options->MaxRays);
621 }
622
623 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
624 int r)
625 {
626 int nvar = P->Dimension - exist - nparam;
627 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
628 return -1;
629 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
630 if (i == -1)
631 return -1;
632 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
633 return -1;
634 return i;
635 }
636
637 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
638 unsigned exist, unsigned nparam, barvinok_options *options)
639 {
640 #ifdef DEBUG_ER
641 fprintf(stderr, "\nER: RedundantRay\n");
642 #endif /* DEBUG_ER */
643
644 Value one;
645 value_init(one);
646 value_set_si(one, 1);
647 int len = P->NbRays-1;
648 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
649 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
650 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
651 for (int j = 0; j < P->NbRays; ++j) {
652 if (j == r)
653 continue;
654 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
655 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
656 }
657
658 P = Rays2Polyhedron(M, options->MaxRays);
659 Matrix_Free(M);
660 evalue *EP = barvinok_enumerate_e_with_options(P, exist, nparam, options);
661 Polyhedron_Free(P);
662 value_clear(one);
663
664 return EP;
665 }
666
667 static evalue* enumerate_redundant_ray(Polyhedron *P,
668 unsigned exist, unsigned nparam, barvinok_options *options)
669 {
670 assert(P->NbBid == 0);
671 int nvar = P->Dimension - exist - nparam;
672 Value m;
673 value_init(m);
674
675 for (int r = 0; r < P->NbRays; ++r) {
676 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
677 continue;
678 int i1 = single_param_pos(P, exist, nparam, r);
679 if (i1 == -1)
680 continue;
681 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
682 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
683 continue;
684 int i2 = single_param_pos(P, exist, nparam, r2);
685 if (i2 == -1)
686 continue;
687 if (i1 != i2)
688 continue;
689
690 value_division(m, P->Ray[r][1+nvar+exist+i1],
691 P->Ray[r2][1+nvar+exist+i1]);
692 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
693 /* r2 divides r => r redundant */
694 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
695 value_clear(m);
696 return enumerate_remove_ray(P, r, exist, nparam, options);
697 }
698
699 value_division(m, P->Ray[r2][1+nvar+exist+i1],
700 P->Ray[r][1+nvar+exist+i1]);
701 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
702 /* r divides r2 => r2 redundant */
703 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
704 value_clear(m);
705 return enumerate_remove_ray(P, r2, exist, nparam, options);
706 }
707 }
708 }
709 value_clear(m);
710 return 0;
711 }
712
713 static Polyhedron *upper_bound(Polyhedron *P,
714 int pos, Value *max, Polyhedron **R)
715 {
716 Value v;
717 int r;
718 value_init(v);
719
720 *R = 0;
721 Polyhedron *N;
722 Polyhedron *B = 0;
723 for (Polyhedron *Q = P; Q; Q = N) {
724 N = Q->next;
725 for (r = 0; r < P->NbRays; ++r) {
726 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
727 value_pos_p(P->Ray[r][1+pos]))
728 break;
729 }
730 if (r < P->NbRays) {
731 Q->next = *R;
732 *R = Q;
733 continue;
734 } else {
735 Q->next = B;
736 B = Q;
737 }
738 for (r = 0; r < P->NbRays; ++r) {
739 if (value_zero_p(P->Ray[r][P->Dimension+1]))
740 continue;
741 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
742 if ((!Q->next && r == 0) || value_gt(v, *max))
743 value_assign(*max, v);
744 }
745 }
746 value_clear(v);
747 return B;
748 }
749
750 static evalue* enumerate_ray(Polyhedron *P,
751 unsigned exist, unsigned nparam, barvinok_options *options)
752 {
753 assert(P->NbBid == 0);
754 int nvar = P->Dimension - exist - nparam;
755
756 int r;
757 for (r = 0; r < P->NbRays; ++r)
758 if (value_zero_p(P->Ray[r][P->Dimension+1]))
759 break;
760 if (r >= P->NbRays)
761 return 0;
762
763 int r2;
764 for (r2 = r+1; r2 < P->NbRays; ++r2)
765 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
766 break;
767 if (r2 < P->NbRays) {
768 if (nvar > 0)
769 return enumerate_sum(P, exist, nparam, options);
770 }
771
772 #ifdef DEBUG_ER
773 fprintf(stderr, "\nER: Ray\n");
774 #endif /* DEBUG_ER */
775
776 Value m;
777 Value one;
778 value_init(m);
779 value_init(one);
780 value_set_si(one, 1);
781 int i = single_param_pos(P, exist, nparam, r);
782 assert(i != -1); // for now;
783
784 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
785 for (int j = 0; j < P->NbRays; ++j) {
786 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
787 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
788 }
789 Polyhedron *S = Rays2Polyhedron(M, options->MaxRays);
790 Matrix_Free(M);
791 Polyhedron *D = DomainDifference(P, S, options->MaxRays);
792 Polyhedron_Free(S);
793 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
794 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
795 Polyhedron *R;
796 D = upper_bound(D, nvar+exist+i, &m, &R);
797 assert(D);
798 Domain_Free(D);
799
800 M = Matrix_Alloc(2, P->Dimension+2);
801 value_set_si(M->p[0][0], 1);
802 value_set_si(M->p[1][0], 1);
803 value_set_si(M->p[0][1+nvar+exist+i], -1);
804 value_set_si(M->p[1][1+nvar+exist+i], 1);
805 value_assign(M->p[0][1+P->Dimension], m);
806 value_oppose(M->p[1][1+P->Dimension], m);
807 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
808 P->Ray[r][1+nvar+exist+i]);
809 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
810 // Matrix_Print(stderr, P_VALUE_FMT, M);
811 D = AddConstraints(M->p[0], 2, P, options->MaxRays);
812 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
813 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
814 P->Ray[r][1+nvar+exist+i]);
815 // Matrix_Print(stderr, P_VALUE_FMT, M);
816 S = AddConstraints(M->p[0], 1, P, options->MaxRays);
817 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
818 Matrix_Free(M);
819
820 evalue *EP = barvinok_enumerate_e_with_options(D, exist, nparam, options);
821 Polyhedron_Free(D);
822 value_clear(one);
823 value_clear(m);
824
825 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
826 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, options->MaxRays);
827 else {
828 M = Matrix_Alloc(1, nparam+2);
829 value_set_si(M->p[0][0], 1);
830 value_set_si(M->p[0][1+i], 1);
831 enumerate_vd_add_ray(EP, M, options->MaxRays);
832 Matrix_Free(M);
833 }
834
835 if (!emptyQ(S)) {
836 evalue *E = barvinok_enumerate_e_with_options(S, exist, nparam, options);
837 eadd(E, EP);
838 evalue_free(E);
839 }
840 Polyhedron_Free(S);
841
842 if (R) {
843 assert(nvar == 0);
844 evalue *ER = enumerate_or(R, exist, nparam, options);
845 eor(ER, EP);
846 free_evalue_refs(ER);
847 free(ER);
848 }
849
850 return EP;
851 }
852
853 static evalue* enumerate_vd(Polyhedron **PA,
854 unsigned exist, unsigned nparam, barvinok_options *options)
855 {
856 Polyhedron *P = *PA;
857 int nvar = P->Dimension - exist - nparam;
858 Param_Polyhedron *PP = NULL;
859 Polyhedron *C = Universe_Polyhedron(nparam);
860 Polyhedron *PR = P;
861 PP = Polyhedron2Param_Polyhedron(PR, C, options);
862 Polyhedron_Free(C);
863
864 int nd;
865 Param_Domain *D, *last;
866 Value c;
867 value_init(c);
868 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
869 ;
870
871 Polyhedron **VD = new Polyhedron *[nd];
872 Polyhedron *TC = true_context(P, C, options->MaxRays);
873 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
874 VD[nd++] = rVD;
875 last = D;
876 END_FORALL_REDUCED_DOMAIN
877 Polyhedron_Free(TC);
878
879 evalue *EP = 0;
880
881 if (nd == 0)
882 EP = evalue_zero();
883
884 /* This doesn't seem to have any effect */
885 if (nd == 1) {
886 Polyhedron *CA = align_context(VD[0], P->Dimension, options->MaxRays);
887 Polyhedron *O = P;
888 P = DomainIntersection(P, CA, options->MaxRays);
889 if (O != *PA)
890 Polyhedron_Free(O);
891 Polyhedron_Free(CA);
892 if (emptyQ(P))
893 EP = evalue_zero();
894 }
895
896 if (PR != *PA)
897 Polyhedron_Free(PR);
898 PR = 0;
899
900 if (!EP && nd > 1) {
901 #ifdef DEBUG_ER
902 fprintf(stderr, "\nER: VD\n");
903 #endif /* DEBUG_ER */
904 for (int i = 0; i < nd; ++i) {
905 Polyhedron *CA = align_context(VD[i], P->Dimension, options->MaxRays);
906 Polyhedron *I = DomainIntersection(P, CA, options->MaxRays);
907
908 if (i == 0)
909 EP = barvinok_enumerate_e_with_options(I, exist, nparam, options);
910 else {
911 evalue *E = barvinok_enumerate_e_with_options(I, exist, nparam,
912 options);
913 eadd(E, EP);
914 evalue_free(E);
915 }
916 Polyhedron_Free(I);
917 Polyhedron_Free(CA);
918 }
919 }
920
921 for (int i = 0; i < nd; ++i)
922 Polyhedron_Free(VD[i]);
923 delete [] VD;
924 value_clear(c);
925
926 if (!EP && nvar == 0) {
927 Value f;
928 value_init(f);
929 Param_Vertices *V, *V2;
930 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
931
932 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
933 bool found = false;
934 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
935 if (V == V2) {
936 found = true;
937 continue;
938 }
939 if (!found)
940 continue;
941 for (int i = 0; i < exist; ++i) {
942 value_oppose(f, V->Vertex->p[i][nparam+1]);
943 Vector_Combine(V->Vertex->p[i],
944 V2->Vertex->p[i],
945 M->p[0] + 1 + nvar + exist,
946 V2->Vertex->p[i][nparam+1],
947 f,
948 nparam+1);
949 int j;
950 for (j = 0; j < nparam; ++j)
951 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
952 break;
953 if (j >= nparam)
954 continue;
955 ConstraintSimplify(M->p[0], M->p[0],
956 P->Dimension+2, &f);
957 value_set_si(M->p[0][0], 0);
958 Polyhedron *para = AddConstraints(M->p[0], 1, P,
959 options->MaxRays);
960 POL_ENSURE_VERTICES(para);
961 if (emptyQ(para)) {
962 Polyhedron_Free(para);
963 continue;
964 }
965 Polyhedron *pos, *neg;
966 value_set_si(M->p[0][0], 1);
967 value_decrement(M->p[0][P->Dimension+1],
968 M->p[0][P->Dimension+1]);
969 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
970 value_set_si(f, -1);
971 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
972 P->Dimension+1);
973 value_decrement(M->p[0][P->Dimension+1],
974 M->p[0][P->Dimension+1]);
975 value_decrement(M->p[0][P->Dimension+1],
976 M->p[0][P->Dimension+1]);
977 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
978 POL_ENSURE_VERTICES(neg);
979 POL_ENSURE_VERTICES(pos);
980 if (emptyQ(neg) && emptyQ(pos)) {
981 Polyhedron_Free(para);
982 Polyhedron_Free(pos);
983 Polyhedron_Free(neg);
984 continue;
985 }
986 #ifdef DEBUG_ER
987 fprintf(stderr, "\nER: Order\n");
988 #endif /* DEBUG_ER */
989 EP = barvinok_enumerate_e_with_options(para, exist, nparam,
990 options);
991 evalue *E;
992 if (!emptyQ(pos)) {
993 E = barvinok_enumerate_e_with_options(pos, exist, nparam,
994 options);
995 eadd(E, EP);
996 evalue_free(E);
997 }
998 if (!emptyQ(neg)) {
999 E = barvinok_enumerate_e_with_options(neg, exist, nparam,
1000 options);
1001 eadd(E, EP);
1002 evalue_free(E);
1003 }
1004 Polyhedron_Free(para);
1005 Polyhedron_Free(pos);
1006 Polyhedron_Free(neg);
1007 break;
1008 }
1009 if (EP)
1010 break;
1011 } END_FORALL_PVertex_in_ParamPolyhedron;
1012 if (EP)
1013 break;
1014 } END_FORALL_PVertex_in_ParamPolyhedron;
1015
1016 if (!EP) {
1017 /* Search for vertex coordinate to split on */
1018 /* First look for one independent of the parameters */
1019 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
1020 for (int i = 0; i < exist; ++i) {
1021 int j;
1022 for (j = 0; j < nparam; ++j)
1023 if (value_notzero_p(V->Vertex->p[i][j]))
1024 break;
1025 if (j < nparam)
1026 continue;
1027 value_set_si(M->p[0][0], 1);
1028 Vector_Set(M->p[0]+1, 0, nvar+exist);
1029 Vector_Copy(V->Vertex->p[i],
1030 M->p[0] + 1 + nvar + exist, nparam+1);
1031 value_oppose(M->p[0][1+nvar+i],
1032 V->Vertex->p[i][nparam+1]);
1033
1034 Polyhedron *pos, *neg;
1035 value_set_si(M->p[0][0], 1);
1036 value_decrement(M->p[0][P->Dimension+1],
1037 M->p[0][P->Dimension+1]);
1038 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
1039 value_set_si(f, -1);
1040 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
1041 P->Dimension+1);
1042 value_decrement(M->p[0][P->Dimension+1],
1043 M->p[0][P->Dimension+1]);
1044 value_decrement(M->p[0][P->Dimension+1],
1045 M->p[0][P->Dimension+1]);
1046 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1047 POL_ENSURE_VERTICES(neg);
1048 POL_ENSURE_VERTICES(pos);
1049 if (emptyQ(neg) || emptyQ(pos)) {
1050 Polyhedron_Free(pos);
1051 Polyhedron_Free(neg);
1052 continue;
1053 }
1054 Polyhedron_Free(pos);
1055 value_increment(M->p[0][P->Dimension+1],
1056 M->p[0][P->Dimension+1]);
1057 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1058 #ifdef DEBUG_ER
1059 fprintf(stderr, "\nER: Vertex\n");
1060 #endif /* DEBUG_ER */
1061 pos->next = neg;
1062 EP = enumerate_or(pos, exist, nparam, options);
1063 break;
1064 }
1065 if (EP)
1066 break;
1067 } END_FORALL_PVertex_in_ParamPolyhedron;
1068 }
1069
1070 if (!EP) {
1071 /* Search for vertex coordinate to split on */
1072 /* Now look for one that depends on the parameters */
1073 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
1074 for (int i = 0; i < exist; ++i) {
1075 value_set_si(M->p[0][0], 1);
1076 Vector_Set(M->p[0]+1, 0, nvar+exist);
1077 Vector_Copy(V->Vertex->p[i],
1078 M->p[0] + 1 + nvar + exist, nparam+1);
1079 value_oppose(M->p[0][1+nvar+i],
1080 V->Vertex->p[i][nparam+1]);
1081
1082 Polyhedron *pos, *neg;
1083 value_set_si(M->p[0][0], 1);
1084 value_decrement(M->p[0][P->Dimension+1],
1085 M->p[0][P->Dimension+1]);
1086 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
1087 value_set_si(f, -1);
1088 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
1089 P->Dimension+1);
1090 value_decrement(M->p[0][P->Dimension+1],
1091 M->p[0][P->Dimension+1]);
1092 value_decrement(M->p[0][P->Dimension+1],
1093 M->p[0][P->Dimension+1]);
1094 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1095 POL_ENSURE_VERTICES(neg);
1096 POL_ENSURE_VERTICES(pos);
1097 if (emptyQ(neg) || emptyQ(pos)) {
1098 Polyhedron_Free(pos);
1099 Polyhedron_Free(neg);
1100 continue;
1101 }
1102 Polyhedron_Free(pos);
1103 value_increment(M->p[0][P->Dimension+1],
1104 M->p[0][P->Dimension+1]);
1105 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1106 #ifdef DEBUG_ER
1107 fprintf(stderr, "\nER: ParamVertex\n");
1108 #endif /* DEBUG_ER */
1109 pos->next = neg;
1110 EP = enumerate_or(pos, exist, nparam, options);
1111 break;
1112 }
1113 if (EP)
1114 break;
1115 } END_FORALL_PVertex_in_ParamPolyhedron;
1116 }
1117
1118 Matrix_Free(M);
1119 value_clear(f);
1120 }
1121
1122 if (PP)
1123 Param_Polyhedron_Free(PP);
1124 *PA = P;
1125
1126 return EP;
1127 }
1128
1129 evalue* barvinok_enumerate_pip(Polyhedron *P, unsigned exist, unsigned nparam,
1130 unsigned MaxRays)
1131 {
1132 evalue *E;
1133 barvinok_options *options = barvinok_options_new_with_defaults();
1134 options->MaxRays = MaxRays;
1135 E = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1136 barvinok_options_free(options);
1137 return E;
1138 }
1139
1140 evalue *barvinok_enumerate_pip_with_options(Polyhedron *P,
1141 unsigned exist, unsigned nparam, struct barvinok_options *options)
1142 {
1143 int nvar = P->Dimension - exist - nparam;
1144 evalue *EP = evalue_zero();
1145 Polyhedron *Q, *N;
1146
1147 #ifdef DEBUG_ER
1148 fprintf(stderr, "\nER: PIP\n");
1149 #endif /* DEBUG_ER */
1150
1151 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
1152 for (Q = D; Q; Q = N) {
1153 N = Q->next;
1154 Q->next = 0;
1155 evalue *E;
1156 exist = Q->Dimension - nvar - nparam;
1157 E = barvinok_enumerate_e_with_options(Q, exist, nparam, options);
1158 Polyhedron_Free(Q);
1159 eadd(E, EP);
1160 evalue_free(E);
1161 }
1162
1163 return EP;
1164 }
1165
1166 evalue *barvinok_enumerate_isl(Polyhedron *P,
1167 unsigned exist, unsigned nparam, struct barvinok_options *options)
1168 {
1169 isl_ctx *ctx = isl_ctx_alloc();
1170 isl_dim *dims;
1171 isl_basic_set *bset;
1172 isl_set *set;
1173 evalue *EP = evalue_zero();
1174 Polyhedron *D, *Q, *N;
1175 Polyhedron *U = Universe_Polyhedron(nparam);
1176
1177 dims = isl_dim_set_alloc(ctx, nparam, P->Dimension - nparam - exist);
1178 bset = isl_basic_set_new_from_polylib(P, dims);
1179
1180 set = isl_basic_set_compute_divs(bset);
1181
1182 D = isl_set_to_polylib(set);
1183 for (Q = D; Q; Q = N) {
1184 N = Q->next;
1185 Q->next = 0;
1186 evalue *E;
1187 E = barvinok_enumerate_with_options(Q, U, options);
1188 Polyhedron_Free(Q);
1189 eadd(E, EP);
1190 evalue_free(E);
1191 }
1192
1193 Polyhedron_Free(U);
1194 isl_set_free(set);
1195 isl_ctx_free(ctx);
1196
1197 return EP;
1198 }
1199
1200 static bool is_single(Value *row, int pos, int len)
1201 {
1202 return First_Non_Zero(row, pos) == -1 &&
1203 First_Non_Zero(row+pos+1, len-pos-1) == -1;
1204 }
1205
1206 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1207 unsigned exist, unsigned nparam, barvinok_options *options);
1208
1209 #ifdef DEBUG_ER
1210 static int er_level = 0;
1211
1212 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1213 unsigned exist, unsigned nparam, barvinok_options *options)
1214 {
1215 fprintf(stderr, "\nER: level %i\n", er_level);
1216
1217 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
1218 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
1219 ++er_level;
1220 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1221 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1222 Polyhedron_Free(P);
1223 --er_level;
1224 return EP;
1225 }
1226 #else
1227 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1228 unsigned exist, unsigned nparam, barvinok_options *options)
1229 {
1230 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1231 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1232 Polyhedron_Free(P);
1233 return EP;
1234 }
1235 #endif
1236
1237 evalue* barvinok_enumerate_e(Polyhedron *P, unsigned exist, unsigned nparam,
1238 unsigned MaxRays)
1239 {
1240 evalue *E;
1241 barvinok_options *options = barvinok_options_new_with_defaults();
1242 options->MaxRays = MaxRays;
1243 E = barvinok_enumerate_e_with_options(P, exist, nparam, options);
1244 barvinok_options_free(options);
1245 return E;
1246 }
1247
1248 static evalue *universal_zero(unsigned nparam)
1249 {
1250 evalue *eres;
1251
1252 eres = ALLOC(evalue);
1253 value_init(eres->d);
1254 value_set_si(eres->d, 0);
1255 eres->x.p = new_enode(partition, 2, nparam);
1256 EVALUE_SET_DOMAIN(eres->x.p->arr[0], Universe_Polyhedron(nparam));
1257 value_set_si(eres->x.p->arr[1].d, 1);
1258 value_init(eres->x.p->arr[1].x.n);
1259
1260 return eres;
1261 }
1262
1263 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1264 unsigned exist, unsigned nparam, barvinok_options *options)
1265 {
1266 if (exist == 0) {
1267 Polyhedron *U = Universe_Polyhedron(nparam);
1268 evalue *EP = barvinok_enumerate_with_options(P, U, options);
1269 //char *param_name[] = {"P", "Q", "R", "S", "T" };
1270 //print_evalue(stdout, EP, param_name);
1271 Polyhedron_Free(U);
1272 return EP;
1273 }
1274
1275 int nvar = P->Dimension - exist - nparam;
1276 int len = P->Dimension + 2;
1277
1278 /* for now */
1279 POL_ENSURE_FACETS(P);
1280 POL_ENSURE_VERTICES(P);
1281
1282 if (emptyQ(P))
1283 return evalue_zero();
1284
1285 if (nvar == 0 && nparam == 0) {
1286 evalue *EP = universal_zero(nparam);
1287 barvinok_count_with_options(P, &EP->x.p->arr[1].x.n, options);
1288 if (value_pos_p(EP->x.p->arr[1].x.n))
1289 value_set_si(EP->x.p->arr[1].x.n, 1);
1290 return EP;
1291 }
1292
1293 int r;
1294 for (r = 0; r < P->NbRays; ++r)
1295 if (value_zero_p(P->Ray[r][0]) ||
1296 value_zero_p(P->Ray[r][P->Dimension+1])) {
1297 int i;
1298 for (i = 0; i < nvar; ++i)
1299 if (value_notzero_p(P->Ray[r][i+1]))
1300 break;
1301 if (i >= nvar)
1302 continue;
1303 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
1304 if (value_notzero_p(P->Ray[r][i+1]))
1305 break;
1306 if (i >= nvar + exist + nparam)
1307 break;
1308 }
1309 if (r < P->NbRays) {
1310 evalue *EP = universal_zero(nparam);
1311 value_set_si(EP->x.p->arr[1].x.n, -1);
1312 return EP;
1313 }
1314
1315 int first;
1316 for (r = 0; r < P->NbEq; ++r)
1317 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
1318 break;
1319 if (r < P->NbEq) {
1320 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
1321 exist-first-1) != -1) {
1322 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1323 #ifdef DEBUG_ER
1324 fprintf(stderr, "\nER: Equality\n");
1325 #endif /* DEBUG_ER */
1326 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1327 options);
1328 Polyhedron_Free(T);
1329 return EP;
1330 } else {
1331 #ifdef DEBUG_ER
1332 fprintf(stderr, "\nER: Fixed\n");
1333 #endif /* DEBUG_ER */
1334 if (first == 0)
1335 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1336 options);
1337 else {
1338 Polyhedron *T = Polyhedron_Copy(P);
1339 Polyhedron_ExchangeColumns(T, nvar+1, nvar+1+first);
1340 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1341 options);
1342 Polyhedron_Free(T);
1343 return EP;
1344 }
1345 }
1346 }
1347
1348 Vector *row = Vector_Alloc(len);
1349 value_set_si(row->p[0], 1);
1350
1351 Value f;
1352 value_init(f);
1353
1354 enum constraint* info = new constraint[exist];
1355 for (int i = 0; i < exist; ++i) {
1356 info[i] = ALL_POS;
1357 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1358 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1359 continue;
1360 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
1361 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1362 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1363 continue;
1364 bool lu_parallel = l_parallel ||
1365 is_single(P->Constraint[u]+nvar+1, i, exist);
1366 value_oppose(f, P->Constraint[u][nvar+i+1]);
1367 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
1368 f, P->Constraint[l][nvar+i+1], len-1);
1369 if (!(info[i] & INDEPENDENT)) {
1370 int j;
1371 for (j = 0; j < exist; ++j)
1372 if (j != i && value_notzero_p(row->p[nvar+j+1]))
1373 break;
1374 if (j == exist) {
1375 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
1376 info[i] = (constraint)(info[i] | INDEPENDENT);
1377 }
1378 }
1379 if (info[i] & ALL_POS) {
1380 value_addto(row->p[len-1], row->p[len-1],
1381 P->Constraint[l][nvar+i+1]);
1382 value_addto(row->p[len-1], row->p[len-1], f);
1383 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1384 value_subtract(row->p[len-1], row->p[len-1], f);
1385 value_decrement(row->p[len-1], row->p[len-1]);
1386 ConstraintSimplify(row->p, row->p, len, &f);
1387 value_set_si(f, -1);
1388 Vector_Scale(row->p+1, row->p+1, f, len-1);
1389 value_decrement(row->p[len-1], row->p[len-1]);
1390 Polyhedron *T = AddConstraints(row->p, 1, P, options->MaxRays);
1391 POL_ENSURE_VERTICES(T);
1392 if (!emptyQ(T)) {
1393 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
1394 info[i] = (constraint)(info[i] ^ ALL_POS);
1395 }
1396 //puts("pos remainder");
1397 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1398 Polyhedron_Free(T);
1399 }
1400 if (!(info[i] & ONE_NEG)) {
1401 if (lu_parallel) {
1402 negative_test_constraint(P->Constraint[l],
1403 P->Constraint[u],
1404 row->p, nvar+i, len, &f);
1405 oppose_constraint(row->p, len, &f);
1406 Polyhedron *T = AddConstraints(row->p, 1, P,
1407 options->MaxRays);
1408 POL_ENSURE_VERTICES(T);
1409 if (emptyQ(T)) {
1410 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
1411 info[i] = (constraint)(info[i] | ONE_NEG);
1412 }
1413 //puts("neg remainder");
1414 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1415 Polyhedron_Free(T);
1416 } else if (!(info[i] & ROT_NEG)) {
1417 if (parallel_constraints(P->Constraint[l],
1418 P->Constraint[u],
1419 row->p, nvar, exist)) {
1420 negative_test_constraint7(P->Constraint[l],
1421 P->Constraint[u],
1422 row->p, nvar, exist,
1423 len, &f);
1424 oppose_constraint(row->p, len, &f);
1425 Polyhedron *T = AddConstraints(row->p, 1, P,
1426 options->MaxRays);
1427 POL_ENSURE_VERTICES(T);
1428 if (emptyQ(T)) {
1429 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
1430 info[i] = (constraint)(info[i] | ROT_NEG);
1431 r = l;
1432 }
1433 //puts("neg remainder");
1434 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1435 Polyhedron_Free(T);
1436 }
1437 }
1438 }
1439 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
1440 goto next;
1441 }
1442 }
1443 if (info[i] & ALL_POS)
1444 break;
1445 next:
1446 ;
1447 }
1448
1449 /*
1450 for (int i = 0; i < exist; ++i)
1451 printf("%i: %i\n", i, info[i]);
1452 */
1453 for (int i = 0; i < exist; ++i)
1454 if (info[i] & ALL_POS) {
1455 #ifdef DEBUG_ER
1456 fprintf(stderr, "\nER: Positive\n");
1457 #endif /* DEBUG_ER */
1458 // Eliminate
1459 // Maybe we should chew off some of the fat here
1460 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
1461 for (int j = 0; j < P->Dimension; ++j)
1462 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
1463 Polyhedron *T = Polyhedron_Image(P, M, options->MaxRays);
1464 Matrix_Free(M);
1465 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1466 options);
1467 Polyhedron_Free(T);
1468 value_clear(f);
1469 Vector_Free(row);
1470 delete [] info;
1471 return EP;
1472 }
1473 for (int i = 0; i < exist; ++i)
1474 if (info[i] & ONE_NEG) {
1475 #ifdef DEBUG_ER
1476 fprintf(stderr, "\nER: Negative\n");
1477 #endif /* DEBUG_ER */
1478 Vector_Free(row);
1479 value_clear(f);
1480 delete [] info;
1481 if (i == 0)
1482 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1483 options);
1484 else {
1485 Polyhedron *T = Polyhedron_Copy(P);
1486 Polyhedron_ExchangeColumns(T, nvar+1, nvar+1+i);
1487 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1488 options);
1489 Polyhedron_Free(T);
1490 return EP;
1491 }
1492 }
1493 for (int i = 0; i < exist; ++i)
1494 if (info[i] & ROT_NEG) {
1495 #ifdef DEBUG_ER
1496 fprintf(stderr, "\nER: Rotate\n");
1497 #endif /* DEBUG_ER */
1498 Vector_Free(row);
1499 value_clear(f);
1500 delete [] info;
1501 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1502 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1503 options);
1504 Polyhedron_Free(T);
1505 return EP;
1506 }
1507 for (int i = 0; i < exist; ++i)
1508 if (info[i] & INDEPENDENT) {
1509 Polyhedron *pos, *neg;
1510
1511 /* Find constraint again and split off negative part */
1512
1513 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1514 row, f, true, &pos, &neg)) {
1515 #ifdef DEBUG_ER
1516 fprintf(stderr, "\nER: Split\n");
1517 #endif /* DEBUG_ER */
1518
1519 evalue *EP =
1520 barvinok_enumerate_e_with_options(neg, exist-1, nparam, options);
1521 evalue *E =
1522 barvinok_enumerate_e_with_options(pos, exist, nparam, options);
1523 eadd(E, EP);
1524 evalue_free(E);
1525 Polyhedron_Free(neg);
1526 Polyhedron_Free(pos);
1527 value_clear(f);
1528 Vector_Free(row);
1529 delete [] info;
1530 return EP;
1531 }
1532 }
1533 delete [] info;
1534
1535 Polyhedron *O = P;
1536 Polyhedron *F;
1537
1538 evalue *EP;
1539
1540 EP = enumerate_line(P, exist, nparam, options);
1541 if (EP)
1542 goto out;
1543
1544 EP = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1545 if (EP)
1546 goto out;
1547
1548 EP = enumerate_redundant_ray(P, exist, nparam, options);
1549 if (EP)
1550 goto out;
1551
1552 EP = enumerate_sure(P, exist, nparam, options);
1553 if (EP)
1554 goto out;
1555
1556 EP = enumerate_ray(P, exist, nparam, options);
1557 if (EP)
1558 goto out;
1559
1560 EP = enumerate_sure2(P, exist, nparam, options);
1561 if (EP)
1562 goto out;
1563
1564 F = unfringe(P, options->MaxRays);
1565 if (!PolyhedronIncludes(F, P)) {
1566 #ifdef DEBUG_ER
1567 fprintf(stderr, "\nER: Fringed\n");
1568 #endif /* DEBUG_ER */
1569 EP = barvinok_enumerate_e_with_options(F, exist, nparam, options);
1570 Polyhedron_Free(F);
1571 goto out;
1572 }
1573 Polyhedron_Free(F);
1574
1575 if (nparam)
1576 EP = enumerate_vd(&P, exist, nparam, options);
1577 if (EP)
1578 goto out2;
1579
1580 if (nvar != 0) {
1581 EP = enumerate_sum(P, exist, nparam, options);
1582 goto out2;
1583 }
1584
1585 assert(nvar == 0);
1586
1587 int i;
1588 Polyhedron *pos, *neg;
1589 for (i = 0; i < exist; ++i)
1590 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1591 row, f, false, &pos, &neg))
1592 break;
1593
1594 assert (i < exist);
1595
1596 pos->next = neg;
1597 EP = enumerate_or(pos, exist, nparam, options);
1598
1599 out2:
1600 if (O != P)
1601 Polyhedron_Free(P);
1602
1603 out:
1604 value_clear(f);
1605 Vector_Free(row);
1606 return EP;
1607 }
lattice point enumeration library