?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue.c: Polyhedron_Insert: add missing return type
[barvinok.git] / barvinok_e.cc
blob35628ec1cbd8c9be35da4373c3c7ee1b7860bfed
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 *CEq;
861 Matrix *CT;
862 Polyhedron *PR = P;
863 PP = Polyhedron2Param_Polyhedron(PR, C, options);
864 Polyhedron_Free(C);
865
866 int nd;
867 Param_Domain *D, *last;
868 Value c;
869 value_init(c);
870 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
871 ;
872
873 Polyhedron **VD = new Polyhedron *[nd];
874 Polyhedron *TC = true_context(P, C, options->MaxRays);
875 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
876 VD[nd++] = rVD;
877 last = D;
878 END_FORALL_REDUCED_DOMAIN
879 Polyhedron_Free(TC);
880
881 evalue *EP = 0;
882
883 if (nd == 0)
884 EP = evalue_zero();
885
886 /* This doesn't seem to have any effect */
887 if (nd == 1) {
888 Polyhedron *CA = align_context(VD[0], P->Dimension, options->MaxRays);
889 Polyhedron *O = P;
890 P = DomainIntersection(P, CA, options->MaxRays);
891 if (O != *PA)
892 Polyhedron_Free(O);
893 Polyhedron_Free(CA);
894 if (emptyQ(P))
895 EP = evalue_zero();
896 }
897
898 if (PR != *PA)
899 Polyhedron_Free(PR);
900 PR = 0;
901
902 if (!EP && nd > 1) {
903 #ifdef DEBUG_ER
904 fprintf(stderr, "\nER: VD\n");
905 #endif /* DEBUG_ER */
906 for (int i = 0; i < nd; ++i) {
907 Polyhedron *CA = align_context(VD[i], P->Dimension, options->MaxRays);
908 Polyhedron *I = DomainIntersection(P, CA, options->MaxRays);
909
910 if (i == 0)
911 EP = barvinok_enumerate_e_with_options(I, exist, nparam, options);
912 else {
913 evalue *E = barvinok_enumerate_e_with_options(I, exist, nparam,
914 options);
915 eadd(E, EP);
916 evalue_free(E);
917 }
918 Polyhedron_Free(I);
919 Polyhedron_Free(CA);
920 }
921 }
922
923 for (int i = 0; i < nd; ++i)
924 Polyhedron_Free(VD[i]);
925 delete [] VD;
926 value_clear(c);
927
928 if (!EP && nvar == 0) {
929 Value f;
930 value_init(f);
931 Param_Vertices *V, *V2;
932 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
933
934 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
935 bool found = false;
936 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
937 if (V == V2) {
938 found = true;
939 continue;
940 }
941 if (!found)
942 continue;
943 for (int i = 0; i < exist; ++i) {
944 value_oppose(f, V->Vertex->p[i][nparam+1]);
945 Vector_Combine(V->Vertex->p[i],
946 V2->Vertex->p[i],
947 M->p[0] + 1 + nvar + exist,
948 V2->Vertex->p[i][nparam+1],
949 f,
950 nparam+1);
951 int j;
952 for (j = 0; j < nparam; ++j)
953 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
954 break;
955 if (j >= nparam)
956 continue;
957 ConstraintSimplify(M->p[0], M->p[0],
958 P->Dimension+2, &f);
959 value_set_si(M->p[0][0], 0);
960 Polyhedron *para = AddConstraints(M->p[0], 1, P,
961 options->MaxRays);
962 POL_ENSURE_VERTICES(para);
963 if (emptyQ(para)) {
964 Polyhedron_Free(para);
965 continue;
966 }
967 Polyhedron *pos, *neg;
968 value_set_si(M->p[0][0], 1);
969 value_decrement(M->p[0][P->Dimension+1],
970 M->p[0][P->Dimension+1]);
971 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
972 value_set_si(f, -1);
973 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
974 P->Dimension+1);
975 value_decrement(M->p[0][P->Dimension+1],
976 M->p[0][P->Dimension+1]);
977 value_decrement(M->p[0][P->Dimension+1],
978 M->p[0][P->Dimension+1]);
979 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
980 POL_ENSURE_VERTICES(neg);
981 POL_ENSURE_VERTICES(pos);
982 if (emptyQ(neg) && emptyQ(pos)) {
983 Polyhedron_Free(para);
984 Polyhedron_Free(pos);
985 Polyhedron_Free(neg);
986 continue;
987 }
988 #ifdef DEBUG_ER
989 fprintf(stderr, "\nER: Order\n");
990 #endif /* DEBUG_ER */
991 EP = barvinok_enumerate_e_with_options(para, exist, nparam,
992 options);
993 evalue *E;
994 if (!emptyQ(pos)) {
995 E = barvinok_enumerate_e_with_options(pos, exist, nparam,
996 options);
997 eadd(E, EP);
998 evalue_free(E);
999 }
1000 if (!emptyQ(neg)) {
1001 E = barvinok_enumerate_e_with_options(neg, exist, nparam,
1002 options);
1003 eadd(E, EP);
1004 evalue_free(E);
1005 }
1006 Polyhedron_Free(para);
1007 Polyhedron_Free(pos);
1008 Polyhedron_Free(neg);
1009 break;
1010 }
1011 if (EP)
1012 break;
1013 } END_FORALL_PVertex_in_ParamPolyhedron;
1014 if (EP)
1015 break;
1016 } END_FORALL_PVertex_in_ParamPolyhedron;
1017
1018 if (!EP) {
1019 /* Search for vertex coordinate to split on */
1020 /* First look for one independent of the parameters */
1021 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
1022 for (int i = 0; i < exist; ++i) {
1023 int j;
1024 for (j = 0; j < nparam; ++j)
1025 if (value_notzero_p(V->Vertex->p[i][j]))
1026 break;
1027 if (j < nparam)
1028 continue;
1029 value_set_si(M->p[0][0], 1);
1030 Vector_Set(M->p[0]+1, 0, nvar+exist);
1031 Vector_Copy(V->Vertex->p[i],
1032 M->p[0] + 1 + nvar + exist, nparam+1);
1033 value_oppose(M->p[0][1+nvar+i],
1034 V->Vertex->p[i][nparam+1]);
1035
1036 Polyhedron *pos, *neg;
1037 value_set_si(M->p[0][0], 1);
1038 value_decrement(M->p[0][P->Dimension+1],
1039 M->p[0][P->Dimension+1]);
1040 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
1041 value_set_si(f, -1);
1042 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
1043 P->Dimension+1);
1044 value_decrement(M->p[0][P->Dimension+1],
1045 M->p[0][P->Dimension+1]);
1046 value_decrement(M->p[0][P->Dimension+1],
1047 M->p[0][P->Dimension+1]);
1048 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1049 POL_ENSURE_VERTICES(neg);
1050 POL_ENSURE_VERTICES(pos);
1051 if (emptyQ(neg) || emptyQ(pos)) {
1052 Polyhedron_Free(pos);
1053 Polyhedron_Free(neg);
1054 continue;
1055 }
1056 Polyhedron_Free(pos);
1057 value_increment(M->p[0][P->Dimension+1],
1058 M->p[0][P->Dimension+1]);
1059 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1060 #ifdef DEBUG_ER
1061 fprintf(stderr, "\nER: Vertex\n");
1062 #endif /* DEBUG_ER */
1063 pos->next = neg;
1064 EP = enumerate_or(pos, exist, nparam, options);
1065 break;
1066 }
1067 if (EP)
1068 break;
1069 } END_FORALL_PVertex_in_ParamPolyhedron;
1070 }
1071
1072 if (!EP) {
1073 /* Search for vertex coordinate to split on */
1074 /* Now look for one that depends on the parameters */
1075 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
1076 for (int i = 0; i < exist; ++i) {
1077 value_set_si(M->p[0][0], 1);
1078 Vector_Set(M->p[0]+1, 0, nvar+exist);
1079 Vector_Copy(V->Vertex->p[i],
1080 M->p[0] + 1 + nvar + exist, nparam+1);
1081 value_oppose(M->p[0][1+nvar+i],
1082 V->Vertex->p[i][nparam+1]);
1083
1084 Polyhedron *pos, *neg;
1085 value_set_si(M->p[0][0], 1);
1086 value_decrement(M->p[0][P->Dimension+1],
1087 M->p[0][P->Dimension+1]);
1088 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
1089 value_set_si(f, -1);
1090 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
1091 P->Dimension+1);
1092 value_decrement(M->p[0][P->Dimension+1],
1093 M->p[0][P->Dimension+1]);
1094 value_decrement(M->p[0][P->Dimension+1],
1095 M->p[0][P->Dimension+1]);
1096 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1097 POL_ENSURE_VERTICES(neg);
1098 POL_ENSURE_VERTICES(pos);
1099 if (emptyQ(neg) || emptyQ(pos)) {
1100 Polyhedron_Free(pos);
1101 Polyhedron_Free(neg);
1102 continue;
1103 }
1104 Polyhedron_Free(pos);
1105 value_increment(M->p[0][P->Dimension+1],
1106 M->p[0][P->Dimension+1]);
1107 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1108 #ifdef DEBUG_ER
1109 fprintf(stderr, "\nER: ParamVertex\n");
1110 #endif /* DEBUG_ER */
1111 pos->next = neg;
1112 EP = enumerate_or(pos, exist, nparam, options);
1113 break;
1114 }
1115 if (EP)
1116 break;
1117 } END_FORALL_PVertex_in_ParamPolyhedron;
1118 }
1119
1120 Matrix_Free(M);
1121 value_clear(f);
1122 }
1123
1124 if (CEq)
1125 Polyhedron_Free(CEq);
1126 if (CT)
1127 Matrix_Free(CT);
1128 if (PP)
1129 Param_Polyhedron_Free(PP);
1130 *PA = P;
1131
1132 return EP;
1133 }
1134
1135 evalue* barvinok_enumerate_pip(Polyhedron *P, unsigned exist, unsigned nparam,
1136 unsigned MaxRays)
1137 {
1138 evalue *E;
1139 barvinok_options *options = barvinok_options_new_with_defaults();
1140 options->MaxRays = MaxRays;
1141 E = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1142 barvinok_options_free(options);
1143 return E;
1144 }
1145
1146 evalue *barvinok_enumerate_pip_with_options(Polyhedron *P,
1147 unsigned exist, unsigned nparam, struct barvinok_options *options)
1148 {
1149 int nvar = P->Dimension - exist - nparam;
1150 evalue *EP = evalue_zero();
1151 Polyhedron *Q, *N;
1152
1153 #ifdef DEBUG_ER
1154 fprintf(stderr, "\nER: PIP\n");
1155 #endif /* DEBUG_ER */
1156
1157 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
1158 for (Q = D; Q; Q = N) {
1159 N = Q->next;
1160 Q->next = 0;
1161 evalue *E;
1162 exist = Q->Dimension - nvar - nparam;
1163 E = barvinok_enumerate_e_with_options(Q, exist, nparam, options);
1164 Polyhedron_Free(Q);
1165 eadd(E, EP);
1166 evalue_free(E);
1167 }
1168
1169 return EP;
1170 }
1171
1172 evalue *barvinok_enumerate_isl(Polyhedron *P,
1173 unsigned exist, unsigned nparam, struct barvinok_options *options)
1174 {
1175 isl_ctx *ctx = isl_ctx_alloc();
1176 isl_dim *dims;
1177 isl_basic_set *bset;
1178 isl_set *set;
1179 evalue *EP = evalue_zero();
1180 Polyhedron *D, *Q, *N;
1181 Polyhedron *U = Universe_Polyhedron(nparam);
1182
1183 dims = isl_dim_set_alloc(ctx, nparam, P->Dimension - nparam - exist);
1184 bset = isl_basic_set_new_from_polylib(P, dims);
1185
1186 set = isl_basic_set_compute_divs(bset);
1187
1188 D = isl_set_to_polylib(set);
1189 for (Q = D; Q; Q = N) {
1190 N = Q->next;
1191 Q->next = 0;
1192 evalue *E;
1193 E = barvinok_enumerate_with_options(Q, U, options);
1194 Polyhedron_Free(Q);
1195 eadd(E, EP);
1196 evalue_free(E);
1197 }
1198
1199 Polyhedron_Free(U);
1200 isl_set_free(set);
1201 isl_ctx_free(ctx);
1202
1203 return EP;
1204 }
1205
1206 static bool is_single(Value *row, int pos, int len)
1207 {
1208 return First_Non_Zero(row, pos) == -1 &&
1209 First_Non_Zero(row+pos+1, len-pos-1) == -1;
1210 }
1211
1212 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1213 unsigned exist, unsigned nparam, barvinok_options *options);
1214
1215 #ifdef DEBUG_ER
1216 static int er_level = 0;
1217
1218 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1219 unsigned exist, unsigned nparam, barvinok_options *options)
1220 {
1221 fprintf(stderr, "\nER: level %i\n", er_level);
1222
1223 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
1224 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
1225 ++er_level;
1226 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1227 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1228 Polyhedron_Free(P);
1229 --er_level;
1230 return EP;
1231 }
1232 #else
1233 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1234 unsigned exist, unsigned nparam, barvinok_options *options)
1235 {
1236 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1237 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1238 Polyhedron_Free(P);
1239 return EP;
1240 }
1241 #endif
1242
1243 evalue* barvinok_enumerate_e(Polyhedron *P, unsigned exist, unsigned nparam,
1244 unsigned MaxRays)
1245 {
1246 evalue *E;
1247 barvinok_options *options = barvinok_options_new_with_defaults();
1248 options->MaxRays = MaxRays;
1249 E = barvinok_enumerate_e_with_options(P, exist, nparam, options);
1250 barvinok_options_free(options);
1251 return E;
1252 }
1253
1254 static evalue *universal_zero(unsigned nparam)
1255 {
1256 evalue *eres;
1257
1258 eres = ALLOC(evalue);
1259 value_init(eres->d);
1260 value_set_si(eres->d, 0);
1261 eres->x.p = new_enode(partition, 2, nparam);
1262 EVALUE_SET_DOMAIN(eres->x.p->arr[0], Universe_Polyhedron(nparam));
1263 value_set_si(eres->x.p->arr[1].d, 1);
1264 value_init(eres->x.p->arr[1].x.n);
1265
1266 return eres;
1267 }
1268
1269 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1270 unsigned exist, unsigned nparam, barvinok_options *options)
1271 {
1272 if (exist == 0) {
1273 Polyhedron *U = Universe_Polyhedron(nparam);
1274 evalue *EP = barvinok_enumerate_with_options(P, U, options);
1275 //char *param_name[] = {"P", "Q", "R", "S", "T" };
1276 //print_evalue(stdout, EP, param_name);
1277 Polyhedron_Free(U);
1278 return EP;
1279 }
1280
1281 int nvar = P->Dimension - exist - nparam;
1282 int len = P->Dimension + 2;
1283
1284 /* for now */
1285 POL_ENSURE_FACETS(P);
1286 POL_ENSURE_VERTICES(P);
1287
1288 if (emptyQ(P))
1289 return evalue_zero();
1290
1291 if (nvar == 0 && nparam == 0) {
1292 evalue *EP = universal_zero(nparam);
1293 barvinok_count_with_options(P, &EP->x.p->arr[1].x.n, options);
1294 if (value_pos_p(EP->x.p->arr[1].x.n))
1295 value_set_si(EP->x.p->arr[1].x.n, 1);
1296 return EP;
1297 }
1298
1299 int r;
1300 for (r = 0; r < P->NbRays; ++r)
1301 if (value_zero_p(P->Ray[r][0]) ||
1302 value_zero_p(P->Ray[r][P->Dimension+1])) {
1303 int i;
1304 for (i = 0; i < nvar; ++i)
1305 if (value_notzero_p(P->Ray[r][i+1]))
1306 break;
1307 if (i >= nvar)
1308 continue;
1309 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
1310 if (value_notzero_p(P->Ray[r][i+1]))
1311 break;
1312 if (i >= nvar + exist + nparam)
1313 break;
1314 }
1315 if (r < P->NbRays) {
1316 evalue *EP = universal_zero(nparam);
1317 value_set_si(EP->x.p->arr[1].x.n, -1);
1318 return EP;
1319 }
1320
1321 int first;
1322 for (r = 0; r < P->NbEq; ++r)
1323 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
1324 break;
1325 if (r < P->NbEq) {
1326 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
1327 exist-first-1) != -1) {
1328 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1329 #ifdef DEBUG_ER
1330 fprintf(stderr, "\nER: Equality\n");
1331 #endif /* DEBUG_ER */
1332 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1333 options);
1334 Polyhedron_Free(T);
1335 return EP;
1336 } else {
1337 #ifdef DEBUG_ER
1338 fprintf(stderr, "\nER: Fixed\n");
1339 #endif /* DEBUG_ER */
1340 if (first == 0)
1341 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1342 options);
1343 else {
1344 Polyhedron *T = Polyhedron_Copy(P);
1345 Polyhedron_ExchangeColumns(T, nvar+1, nvar+1+first);
1346 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1347 options);
1348 Polyhedron_Free(T);
1349 return EP;
1350 }
1351 }
1352 }
1353
1354 Vector *row = Vector_Alloc(len);
1355 value_set_si(row->p[0], 1);
1356
1357 Value f;
1358 value_init(f);
1359
1360 enum constraint* info = new constraint[exist];
1361 for (int i = 0; i < exist; ++i) {
1362 info[i] = ALL_POS;
1363 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1364 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1365 continue;
1366 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
1367 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1368 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1369 continue;
1370 bool lu_parallel = l_parallel ||
1371 is_single(P->Constraint[u]+nvar+1, i, exist);
1372 value_oppose(f, P->Constraint[u][nvar+i+1]);
1373 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
1374 f, P->Constraint[l][nvar+i+1], len-1);
1375 if (!(info[i] & INDEPENDENT)) {
1376 int j;
1377 for (j = 0; j < exist; ++j)
1378 if (j != i && value_notzero_p(row->p[nvar+j+1]))
1379 break;
1380 if (j == exist) {
1381 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
1382 info[i] = (constraint)(info[i] | INDEPENDENT);
1383 }
1384 }
1385 if (info[i] & ALL_POS) {
1386 value_addto(row->p[len-1], row->p[len-1],
1387 P->Constraint[l][nvar+i+1]);
1388 value_addto(row->p[len-1], row->p[len-1], f);
1389 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1390 value_subtract(row->p[len-1], row->p[len-1], f);
1391 value_decrement(row->p[len-1], row->p[len-1]);
1392 ConstraintSimplify(row->p, row->p, len, &f);
1393 value_set_si(f, -1);
1394 Vector_Scale(row->p+1, row->p+1, f, len-1);
1395 value_decrement(row->p[len-1], row->p[len-1]);
1396 Polyhedron *T = AddConstraints(row->p, 1, P, options->MaxRays);
1397 POL_ENSURE_VERTICES(T);
1398 if (!emptyQ(T)) {
1399 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
1400 info[i] = (constraint)(info[i] ^ ALL_POS);
1401 }
1402 //puts("pos remainder");
1403 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1404 Polyhedron_Free(T);
1405 }
1406 if (!(info[i] & ONE_NEG)) {
1407 if (lu_parallel) {
1408 negative_test_constraint(P->Constraint[l],
1409 P->Constraint[u],
1410 row->p, nvar+i, len, &f);
1411 oppose_constraint(row->p, len, &f);
1412 Polyhedron *T = AddConstraints(row->p, 1, P,
1413 options->MaxRays);
1414 POL_ENSURE_VERTICES(T);
1415 if (emptyQ(T)) {
1416 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
1417 info[i] = (constraint)(info[i] | ONE_NEG);
1418 }
1419 //puts("neg remainder");
1420 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1421 Polyhedron_Free(T);
1422 } else if (!(info[i] & ROT_NEG)) {
1423 if (parallel_constraints(P->Constraint[l],
1424 P->Constraint[u],
1425 row->p, nvar, exist)) {
1426 negative_test_constraint7(P->Constraint[l],
1427 P->Constraint[u],
1428 row->p, nvar, exist,
1429 len, &f);
1430 oppose_constraint(row->p, len, &f);
1431 Polyhedron *T = AddConstraints(row->p, 1, P,
1432 options->MaxRays);
1433 POL_ENSURE_VERTICES(T);
1434 if (emptyQ(T)) {
1435 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
1436 info[i] = (constraint)(info[i] | ROT_NEG);
1437 r = l;
1438 }
1439 //puts("neg remainder");
1440 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1441 Polyhedron_Free(T);
1442 }
1443 }
1444 }
1445 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
1446 goto next;
1447 }
1448 }
1449 if (info[i] & ALL_POS)
1450 break;
1451 next:
1452 ;
1453 }
1454
1455 /*
1456 for (int i = 0; i < exist; ++i)
1457 printf("%i: %i\n", i, info[i]);
1458 */
1459 for (int i = 0; i < exist; ++i)
1460 if (info[i] & ALL_POS) {
1461 #ifdef DEBUG_ER
1462 fprintf(stderr, "\nER: Positive\n");
1463 #endif /* DEBUG_ER */
1464 // Eliminate
1465 // Maybe we should chew off some of the fat here
1466 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
1467 for (int j = 0; j < P->Dimension; ++j)
1468 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
1469 Polyhedron *T = Polyhedron_Image(P, M, options->MaxRays);
1470 Matrix_Free(M);
1471 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1472 options);
1473 Polyhedron_Free(T);
1474 value_clear(f);
1475 Vector_Free(row);
1476 delete [] info;
1477 return EP;
1478 }
1479 for (int i = 0; i < exist; ++i)
1480 if (info[i] & ONE_NEG) {
1481 #ifdef DEBUG_ER
1482 fprintf(stderr, "\nER: Negative\n");
1483 #endif /* DEBUG_ER */
1484 Vector_Free(row);
1485 value_clear(f);
1486 delete [] info;
1487 if (i == 0)
1488 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1489 options);
1490 else {
1491 Polyhedron *T = Polyhedron_Copy(P);
1492 Polyhedron_ExchangeColumns(T, nvar+1, nvar+1+i);
1493 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1494 options);
1495 Polyhedron_Free(T);
1496 return EP;
1497 }
1498 }
1499 for (int i = 0; i < exist; ++i)
1500 if (info[i] & ROT_NEG) {
1501 #ifdef DEBUG_ER
1502 fprintf(stderr, "\nER: Rotate\n");
1503 #endif /* DEBUG_ER */
1504 Vector_Free(row);
1505 value_clear(f);
1506 delete [] info;
1507 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1508 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1509 options);
1510 Polyhedron_Free(T);
1511 return EP;
1512 }
1513 for (int i = 0; i < exist; ++i)
1514 if (info[i] & INDEPENDENT) {
1515 Polyhedron *pos, *neg;
1516
1517 /* Find constraint again and split off negative part */
1518
1519 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1520 row, f, true, &pos, &neg)) {
1521 #ifdef DEBUG_ER
1522 fprintf(stderr, "\nER: Split\n");
1523 #endif /* DEBUG_ER */
1524
1525 evalue *EP =
1526 barvinok_enumerate_e_with_options(neg, exist-1, nparam, options);
1527 evalue *E =
1528 barvinok_enumerate_e_with_options(pos, exist, nparam, options);
1529 eadd(E, EP);
1530 evalue_free(E);
1531 Polyhedron_Free(neg);
1532 Polyhedron_Free(pos);
1533 value_clear(f);
1534 Vector_Free(row);
1535 delete [] info;
1536 return EP;
1537 }
1538 }
1539 delete [] info;
1540
1541 Polyhedron *O = P;
1542 Polyhedron *F;
1543
1544 evalue *EP;
1545
1546 EP = enumerate_line(P, exist, nparam, options);
1547 if (EP)
1548 goto out;
1549
1550 EP = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1551 if (EP)
1552 goto out;
1553
1554 EP = enumerate_redundant_ray(P, exist, nparam, options);
1555 if (EP)
1556 goto out;
1557
1558 EP = enumerate_sure(P, exist, nparam, options);
1559 if (EP)
1560 goto out;
1561
1562 EP = enumerate_ray(P, exist, nparam, options);
1563 if (EP)
1564 goto out;
1565
1566 EP = enumerate_sure2(P, exist, nparam, options);
1567 if (EP)
1568 goto out;
1569
1570 F = unfringe(P, options->MaxRays);
1571 if (!PolyhedronIncludes(F, P)) {
1572 #ifdef DEBUG_ER
1573 fprintf(stderr, "\nER: Fringed\n");
1574 #endif /* DEBUG_ER */
1575 EP = barvinok_enumerate_e_with_options(F, exist, nparam, options);
1576 Polyhedron_Free(F);
1577 goto out;
1578 }
1579 Polyhedron_Free(F);
1580
1581 if (nparam)
1582 EP = enumerate_vd(&P, exist, nparam, options);
1583 if (EP)
1584 goto out2;
1585
1586 if (nvar != 0) {
1587 EP = enumerate_sum(P, exist, nparam, options);
1588 goto out2;
1589 }
1590
1591 assert(nvar == 0);
1592
1593 int i;
1594 Polyhedron *pos, *neg;
1595 for (i = 0; i < exist; ++i)
1596 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1597 row, f, false, &pos, &neg))
1598 break;
1599
1600 assert (i < exist);
1601
1602 pos->next = neg;
1603 EP = enumerate_or(pos, exist, nparam, options);
1604
1605 out2:
1606 if (O != P)
1607 Polyhedron_Free(P);
1608
1609 out:
1610 value_clear(f);
1611 Vector_Free(row);
1612 return EP;
1613 }
lattice point enumeration library