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