?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'topcom'
[barvinok.git] / barvinok_e.cc
bloba788f9361a0a8f0e7c9a68e77fb028bf38c22a8e
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 evalue *barvinok_enumerate_pip_with_options(Polyhedron *P,
1163 unsigned exist, unsigned nparam, struct barvinok_options *options)
1164 {
1165 int nvar = P->Dimension - exist - nparam;
1166 evalue *EP = evalue_zero();
1167 Polyhedron *Q, *N;
1168
1169 #ifdef DEBUG_ER
1170 fprintf(stderr, "\nER: PIP\n");
1171 #endif /* DEBUG_ER */
1172
1173 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
1174 for (Q = D; Q; Q = N) {
1175 N = Q->next;
1176 Q->next = 0;
1177 evalue *E;
1178 exist = Q->Dimension - nvar - nparam;
1179 E = barvinok_enumerate_e_with_options(Q, exist, nparam, options);
1180 Polyhedron_Free(Q);
1181 eadd(E, EP);
1182 free_evalue_refs(E);
1183 free(E);
1184 }
1185
1186 return EP;
1187 }
1188
1189 static bool is_single(Value *row, int pos, int len)
1190 {
1191 return First_Non_Zero(row, pos) == -1 &&
1192 First_Non_Zero(row+pos+1, len-pos-1) == -1;
1193 }
1194
1195 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1196 unsigned exist, unsigned nparam, barvinok_options *options);
1197
1198 #ifdef DEBUG_ER
1199 static int er_level = 0;
1200
1201 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1202 unsigned exist, unsigned nparam, barvinok_options *options)
1203 {
1204 fprintf(stderr, "\nER: level %i\n", er_level);
1205
1206 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
1207 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
1208 ++er_level;
1209 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1210 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1211 Polyhedron_Free(P);
1212 --er_level;
1213 return EP;
1214 }
1215 #else
1216 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1217 unsigned exist, unsigned nparam, barvinok_options *options)
1218 {
1219 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1220 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1221 Polyhedron_Free(P);
1222 return EP;
1223 }
1224 #endif
1225
1226 evalue* barvinok_enumerate_e(Polyhedron *P, unsigned exist, unsigned nparam,
1227 unsigned MaxRays)
1228 {
1229 evalue *E;
1230 barvinok_options *options = barvinok_options_new_with_defaults();
1231 options->MaxRays = MaxRays;
1232 E = barvinok_enumerate_e_with_options(P, exist, nparam, options);
1233 barvinok_options_free(options);
1234 return E;
1235 }
1236
1237 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1238 unsigned exist, unsigned nparam, barvinok_options *options)
1239 {
1240 if (exist == 0) {
1241 Polyhedron *U = Universe_Polyhedron(nparam);
1242 evalue *EP = barvinok_enumerate_with_options(P, U, options);
1243 //char *param_name[] = {"P", "Q", "R", "S", "T" };
1244 //print_evalue(stdout, EP, param_name);
1245 Polyhedron_Free(U);
1246 return EP;
1247 }
1248
1249 int nvar = P->Dimension - exist - nparam;
1250 int len = P->Dimension + 2;
1251
1252 /* for now */
1253 POL_ENSURE_FACETS(P);
1254 POL_ENSURE_VERTICES(P);
1255
1256 if (emptyQ(P))
1257 return evalue_zero();
1258
1259 if (nvar == 0 && nparam == 0) {
1260 evalue *EP = evalue_zero();
1261 barvinok_count_with_options(P, &EP->x.n, options);
1262 if (value_pos_p(EP->x.n))
1263 value_set_si(EP->x.n, 1);
1264 return EP;
1265 }
1266
1267 int r;
1268 for (r = 0; r < P->NbRays; ++r)
1269 if (value_zero_p(P->Ray[r][0]) ||
1270 value_zero_p(P->Ray[r][P->Dimension+1])) {
1271 int i;
1272 for (i = 0; i < nvar; ++i)
1273 if (value_notzero_p(P->Ray[r][i+1]))
1274 break;
1275 if (i >= nvar)
1276 continue;
1277 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
1278 if (value_notzero_p(P->Ray[r][i+1]))
1279 break;
1280 if (i >= nvar + exist + nparam)
1281 break;
1282 }
1283 if (r < P->NbRays) {
1284 evalue *EP = evalue_zero();
1285 value_set_si(EP->x.n, -1);
1286 return EP;
1287 }
1288
1289 int first;
1290 for (r = 0; r < P->NbEq; ++r)
1291 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
1292 break;
1293 if (r < P->NbEq) {
1294 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
1295 exist-first-1) != -1) {
1296 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1297 #ifdef DEBUG_ER
1298 fprintf(stderr, "\nER: Equality\n");
1299 #endif /* DEBUG_ER */
1300 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1301 options);
1302 Polyhedron_Free(T);
1303 return EP;
1304 } else {
1305 #ifdef DEBUG_ER
1306 fprintf(stderr, "\nER: Fixed\n");
1307 #endif /* DEBUG_ER */
1308 if (first == 0)
1309 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1310 options);
1311 else {
1312 Polyhedron *T = Polyhedron_Copy(P);
1313 SwapColumns(T, nvar+1, nvar+1+first);
1314 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1315 options);
1316 Polyhedron_Free(T);
1317 return EP;
1318 }
1319 }
1320 }
1321
1322 Vector *row = Vector_Alloc(len);
1323 value_set_si(row->p[0], 1);
1324
1325 Value f;
1326 value_init(f);
1327
1328 enum constraint* info = new constraint[exist];
1329 for (int i = 0; i < exist; ++i) {
1330 info[i] = ALL_POS;
1331 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1332 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1333 continue;
1334 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
1335 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1336 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1337 continue;
1338 bool lu_parallel = l_parallel ||
1339 is_single(P->Constraint[u]+nvar+1, i, exist);
1340 value_oppose(f, P->Constraint[u][nvar+i+1]);
1341 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
1342 f, P->Constraint[l][nvar+i+1], len-1);
1343 if (!(info[i] & INDEPENDENT)) {
1344 int j;
1345 for (j = 0; j < exist; ++j)
1346 if (j != i && value_notzero_p(row->p[nvar+j+1]))
1347 break;
1348 if (j == exist) {
1349 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
1350 info[i] = (constraint)(info[i] | INDEPENDENT);
1351 }
1352 }
1353 if (info[i] & ALL_POS) {
1354 value_addto(row->p[len-1], row->p[len-1],
1355 P->Constraint[l][nvar+i+1]);
1356 value_addto(row->p[len-1], row->p[len-1], f);
1357 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1358 value_subtract(row->p[len-1], row->p[len-1], f);
1359 value_decrement(row->p[len-1], row->p[len-1]);
1360 ConstraintSimplify(row->p, row->p, len, &f);
1361 value_set_si(f, -1);
1362 Vector_Scale(row->p+1, row->p+1, f, len-1);
1363 value_decrement(row->p[len-1], row->p[len-1]);
1364 Polyhedron *T = AddConstraints(row->p, 1, P, options->MaxRays);
1365 POL_ENSURE_VERTICES(T);
1366 if (!emptyQ(T)) {
1367 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
1368 info[i] = (constraint)(info[i] ^ ALL_POS);
1369 }
1370 //puts("pos remainder");
1371 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1372 Polyhedron_Free(T);
1373 }
1374 if (!(info[i] & ONE_NEG)) {
1375 if (lu_parallel) {
1376 negative_test_constraint(P->Constraint[l],
1377 P->Constraint[u],
1378 row->p, nvar+i, len, &f);
1379 oppose_constraint(row->p, len, &f);
1380 Polyhedron *T = AddConstraints(row->p, 1, P,
1381 options->MaxRays);
1382 POL_ENSURE_VERTICES(T);
1383 if (emptyQ(T)) {
1384 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
1385 info[i] = (constraint)(info[i] | ONE_NEG);
1386 }
1387 //puts("neg remainder");
1388 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1389 Polyhedron_Free(T);
1390 } else if (!(info[i] & ROT_NEG)) {
1391 if (parallel_constraints(P->Constraint[l],
1392 P->Constraint[u],
1393 row->p, nvar, exist)) {
1394 negative_test_constraint7(P->Constraint[l],
1395 P->Constraint[u],
1396 row->p, nvar, exist,
1397 len, &f);
1398 oppose_constraint(row->p, len, &f);
1399 Polyhedron *T = AddConstraints(row->p, 1, P,
1400 options->MaxRays);
1401 POL_ENSURE_VERTICES(T);
1402 if (emptyQ(T)) {
1403 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
1404 info[i] = (constraint)(info[i] | ROT_NEG);
1405 r = l;
1406 }
1407 //puts("neg remainder");
1408 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1409 Polyhedron_Free(T);
1410 }
1411 }
1412 }
1413 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
1414 goto next;
1415 }
1416 }
1417 if (info[i] & ALL_POS)
1418 break;
1419 next:
1420 ;
1421 }
1422
1423 /*
1424 for (int i = 0; i < exist; ++i)
1425 printf("%i: %i\n", i, info[i]);
1426 */
1427 for (int i = 0; i < exist; ++i)
1428 if (info[i] & ALL_POS) {
1429 #ifdef DEBUG_ER
1430 fprintf(stderr, "\nER: Positive\n");
1431 #endif /* DEBUG_ER */
1432 // Eliminate
1433 // Maybe we should chew off some of the fat here
1434 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
1435 for (int j = 0; j < P->Dimension; ++j)
1436 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
1437 Polyhedron *T = Polyhedron_Image(P, M, options->MaxRays);
1438 Matrix_Free(M);
1439 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1440 options);
1441 Polyhedron_Free(T);
1442 value_clear(f);
1443 Vector_Free(row);
1444 delete [] info;
1445 return EP;
1446 }
1447 for (int i = 0; i < exist; ++i)
1448 if (info[i] & ONE_NEG) {
1449 #ifdef DEBUG_ER
1450 fprintf(stderr, "\nER: Negative\n");
1451 #endif /* DEBUG_ER */
1452 Vector_Free(row);
1453 value_clear(f);
1454 delete [] info;
1455 if (i == 0)
1456 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1457 options);
1458 else {
1459 Polyhedron *T = Polyhedron_Copy(P);
1460 SwapColumns(T, nvar+1, nvar+1+i);
1461 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1462 options);
1463 Polyhedron_Free(T);
1464 return EP;
1465 }
1466 }
1467 for (int i = 0; i < exist; ++i)
1468 if (info[i] & ROT_NEG) {
1469 #ifdef DEBUG_ER
1470 fprintf(stderr, "\nER: Rotate\n");
1471 #endif /* DEBUG_ER */
1472 Vector_Free(row);
1473 value_clear(f);
1474 delete [] info;
1475 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1476 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1477 options);
1478 Polyhedron_Free(T);
1479 return EP;
1480 }
1481 for (int i = 0; i < exist; ++i)
1482 if (info[i] & INDEPENDENT) {
1483 Polyhedron *pos, *neg;
1484
1485 /* Find constraint again and split off negative part */
1486
1487 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1488 row, f, true, &pos, &neg)) {
1489 #ifdef DEBUG_ER
1490 fprintf(stderr, "\nER: Split\n");
1491 #endif /* DEBUG_ER */
1492
1493 evalue *EP =
1494 barvinok_enumerate_e_with_options(neg, exist-1, nparam, options);
1495 evalue *E =
1496 barvinok_enumerate_e_with_options(pos, exist, nparam, options);
1497 eadd(E, EP);
1498 free_evalue_refs(E);
1499 free(E);
1500 Polyhedron_Free(neg);
1501 Polyhedron_Free(pos);
1502 value_clear(f);
1503 Vector_Free(row);
1504 delete [] info;
1505 return EP;
1506 }
1507 }
1508 delete [] info;
1509
1510 Polyhedron *O = P;
1511 Polyhedron *F;
1512
1513 evalue *EP;
1514
1515 EP = enumerate_line(P, exist, nparam, options);
1516 if (EP)
1517 goto out;
1518
1519 EP = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1520 if (EP)
1521 goto out;
1522
1523 EP = enumerate_redundant_ray(P, exist, nparam, options);
1524 if (EP)
1525 goto out;
1526
1527 EP = enumerate_sure(P, exist, nparam, options);
1528 if (EP)
1529 goto out;
1530
1531 EP = enumerate_ray(P, exist, nparam, options);
1532 if (EP)
1533 goto out;
1534
1535 EP = enumerate_sure2(P, exist, nparam, options);
1536 if (EP)
1537 goto out;
1538
1539 F = unfringe(P, options->MaxRays);
1540 if (!PolyhedronIncludes(F, P)) {
1541 #ifdef DEBUG_ER
1542 fprintf(stderr, "\nER: Fringed\n");
1543 #endif /* DEBUG_ER */
1544 EP = barvinok_enumerate_e_with_options(F, exist, nparam, options);
1545 Polyhedron_Free(F);
1546 goto out;
1547 }
1548 Polyhedron_Free(F);
1549
1550 if (nparam)
1551 EP = enumerate_vd(&P, exist, nparam, options);
1552 if (EP)
1553 goto out2;
1554
1555 if (nvar != 0) {
1556 EP = enumerate_sum(P, exist, nparam, options);
1557 goto out2;
1558 }
1559
1560 assert(nvar == 0);
1561
1562 int i;
1563 Polyhedron *pos, *neg;
1564 for (i = 0; i < exist; ++i)
1565 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1566 row, f, false, &pos, &neg))
1567 break;
1568
1569 assert (i < exist);
1570
1571 pos->next = neg;
1572 EP = enumerate_or(pos, exist, nparam, options);
1573
1574 out2:
1575 if (O != P)
1576 Polyhedron_Free(P);
1577
1578 out:
1579 value_clear(f);
1580 Vector_Free(row);
1581 return EP;
1582 }
lattice point enumeration library