?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue.c: eadd: handle some special cases
[barvinok.git] / barvinok_e.cc
blobc4dd6a3688eaeaaf2514c1bf6ae44d028da630ca
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 evalue_free(EN);
368 Polyhedron_Free(D);
369 }
370
371 reduce_evalue(EP);
372
373 return EP;
374 }
375
376 static evalue* enumerate_sum(Polyhedron *P,
377 unsigned exist, unsigned nparam, barvinok_options *options)
378 {
379 int nvar = P->Dimension - exist - nparam;
380 int toswap = nvar < exist ? nvar : exist;
381 for (int i = 0; i < toswap; ++i)
382 SwapColumns(P, 1 + i, nvar+exist - i);
383 nparam += nvar;
384
385 #ifdef DEBUG_ER
386 fprintf(stderr, "\nER: Sum\n");
387 #endif /* DEBUG_ER */
388
389 evalue *EP = barvinok_enumerate_e_with_options(P, exist, nparam, options);
390
391 evalue_split_domains_into_orthants(EP, options->MaxRays);
392 reduce_evalue(EP);
393 evalue_range_reduction(EP);
394
395 evalue_frac2floor(EP);
396
397 evalue *sum = evalue_sum(EP, nvar, options->MaxRays);
398
399 evalue_free(EP);
400 EP = sum;
401
402 evalue_range_reduction(EP);
403
404 return EP;
405 }
406
407 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
408 unsigned exist, unsigned nparam, barvinok_options *options)
409 {
410 int nvar = P->Dimension - exist - nparam;
411
412 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
413 for (int i = 0; i < exist; ++i)
414 value_set_si(M->p[i][nvar+i+1], 1);
415 Polyhedron *O = S;
416 S = DomainAddRays(S, M, options->MaxRays);
417 Polyhedron_Free(O);
418 Polyhedron *F = DomainAddRays(P, M, options->MaxRays);
419 Polyhedron *D = DomainDifference(F, S, options->MaxRays);
420 O = D;
421 D = Disjoint_Domain(D, 0, options->MaxRays);
422 Polyhedron_Free(F);
423 Domain_Free(O);
424 Matrix_Free(M);
425
426 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
427 for (int j = 0; j < nvar; ++j)
428 value_set_si(M->p[j][j], 1);
429 for (int j = 0; j < nparam+1; ++j)
430 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
431 Polyhedron *T = Polyhedron_Image(S, M, options->MaxRays);
432 evalue *EP = barvinok_enumerate_e_with_options(T, 0, nparam, options);
433 Polyhedron_Free(S);
434 Polyhedron_Free(T);
435 Matrix_Free(M);
436
437 for (Polyhedron *Q = D; Q; Q = Q->next) {
438 Polyhedron *N = Q->next;
439 Q->next = 0;
440 T = DomainIntersection(P, Q, options->MaxRays);
441 evalue *E = barvinok_enumerate_e_with_options(T, exist, nparam, options);
442 eadd(E, EP);
443 evalue_free(E);
444 Polyhedron_Free(T);
445 Q->next = N;
446 }
447 Domain_Free(D);
448 return EP;
449 }
450
451 static evalue* enumerate_sure(Polyhedron *P,
452 unsigned exist, unsigned nparam, barvinok_options *options)
453 {
454 int i;
455 Polyhedron *S = P;
456 int nvar = P->Dimension - exist - nparam;
457 Value lcm;
458 Value f;
459 value_init(lcm);
460 value_init(f);
461
462 for (i = 0; i < exist; ++i) {
463 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
464 int c = 0;
465 value_set_si(lcm, 1);
466 for (int j = 0; j < S->NbConstraints; ++j) {
467 if (value_negz_p(S->Constraint[j][1+nvar+i]))
468 continue;
469 if (value_one_p(S->Constraint[j][1+nvar+i]))
470 continue;
471 value_lcm(lcm, lcm, S->Constraint[j][1+nvar+i]);
472 }
473
474 for (int j = 0; j < S->NbConstraints; ++j) {
475 if (value_negz_p(S->Constraint[j][1+nvar+i]))
476 continue;
477 if (value_one_p(S->Constraint[j][1+nvar+i]))
478 continue;
479 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
480 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
481 value_subtract(M->p[c][S->Dimension+1],
482 M->p[c][S->Dimension+1],
483 lcm);
484 value_increment(M->p[c][S->Dimension+1],
485 M->p[c][S->Dimension+1]);
486 ++c;
487 }
488 Polyhedron *O = S;
489 S = AddConstraints(M->p[0], c, S, options->MaxRays);
490 if (O != P)
491 Polyhedron_Free(O);
492 Matrix_Free(M);
493 if (emptyQ(S)) {
494 Polyhedron_Free(S);
495 value_clear(lcm);
496 value_clear(f);
497 return 0;
498 }
499 }
500 value_clear(lcm);
501 value_clear(f);
502
503 #ifdef DEBUG_ER
504 fprintf(stderr, "\nER: Sure\n");
505 #endif /* DEBUG_ER */
506
507 return split_sure(P, S, exist, nparam, options);
508 }
509
510 static evalue* enumerate_sure2(Polyhedron *P,
511 unsigned exist, unsigned nparam, barvinok_options *options)
512 {
513 int nvar = P->Dimension - exist - nparam;
514 int r;
515 for (r = 0; r < P->NbRays; ++r)
516 if (value_one_p(P->Ray[r][0]) &&
517 value_one_p(P->Ray[r][P->Dimension+1]))
518 break;
519
520 if (r >= P->NbRays)
521 return 0;
522
523 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
524 for (int i = 0; i < nvar; ++i)
525 value_set_si(M->p[i][1+i], 1);
526 for (int i = 0; i < nparam; ++i)
527 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
528 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
529 value_set_si(M->p[nvar+nparam][0], 1);
530 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
531 Polyhedron * F = Rays2Polyhedron(M, options->MaxRays);
532 Matrix_Free(M);
533
534 Polyhedron *I = DomainIntersection(F, P, options->MaxRays);
535 Polyhedron_Free(F);
536
537 #ifdef DEBUG_ER
538 fprintf(stderr, "\nER: Sure2\n");
539 #endif /* DEBUG_ER */
540
541 return split_sure(P, I, exist, nparam, options);
542 }
543
544 static evalue* enumerate_cyclic(Polyhedron *P,
545 unsigned exist, unsigned nparam,
546 evalue * EP, int r, int p, unsigned MaxRays)
547 {
548 int nvar = P->Dimension - exist - nparam;
549
550 /* If EP in its fractional maps only contains references
551 * to the remainder parameter with appropriate coefficients
552 * then we could in principle avoid adding existentially
553 * quantified variables to the validity domains.
554 * We'd have to replace the remainder by m { p/m }
555 * and multiply with an appropriate factor that is one
556 * only in the appropriate range.
557 * This last multiplication can be avoided if EP
558 * has a single validity domain with no (further)
559 * constraints on the remainder parameter
560 */
561
562 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
563 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
564 for (int j = 0; j < nparam; ++j)
565 if (j != p)
566 value_set_si(CT->p[j][j], 1);
567 value_set_si(CT->p[p][nparam+1], 1);
568 value_set_si(CT->p[nparam][nparam+2], 1);
569 value_set_si(M->p[0][1+p], -1);
570 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
571 value_set_si(M->p[0][1+nparam+1], 1);
572 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
573 Matrix_Free(M);
574 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
575 Polyhedron_Free(CEq);
576 Matrix_Free(CT);
577
578 return EP;
579 }
580
581 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
582 {
583 if (value_notzero_p(EP->d))
584 return;
585
586 assert(EP->x.p->type == partition);
587 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
588 for (int i = 0; i < EP->x.p->size/2; ++i) {
589 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
590 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
591 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
592 Domain_Free(D);
593 }
594 }
595
596 static evalue* enumerate_line(Polyhedron *P,
597 unsigned exist, unsigned nparam, barvinok_options *options)
598 {
599 if (P->NbBid == 0)
600 return 0;
601
602 #ifdef DEBUG_ER
603 fprintf(stderr, "\nER: Line\n");
604 #endif /* DEBUG_ER */
605
606 int nvar = P->Dimension - exist - nparam;
607 int i, j;
608 for (i = 0; i < nparam; ++i)
609 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
610 break;
611 assert(i < nparam);
612 for (j = i+1; j < nparam; ++j)
613 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
614 break;
615 assert(j >= nparam); // for now
616
617 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
618 value_set_si(M->p[0][0], 1);
619 value_set_si(M->p[0][1+nvar+exist+i], 1);
620 value_set_si(M->p[1][0], 1);
621 value_set_si(M->p[1][1+nvar+exist+i], -1);
622 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
623 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
624 Polyhedron *S = AddConstraints(M->p[0], 2, P, options->MaxRays);
625 evalue *EP = barvinok_enumerate_e_with_options(S, exist, nparam, options);
626 Polyhedron_Free(S);
627 Matrix_Free(M);
628
629 return enumerate_cyclic(P, exist, nparam, EP, 0, i, options->MaxRays);
630 }
631
632 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
633 int r)
634 {
635 int nvar = P->Dimension - exist - nparam;
636 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
637 return -1;
638 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
639 if (i == -1)
640 return -1;
641 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
642 return -1;
643 return i;
644 }
645
646 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
647 unsigned exist, unsigned nparam, barvinok_options *options)
648 {
649 #ifdef DEBUG_ER
650 fprintf(stderr, "\nER: RedundantRay\n");
651 #endif /* DEBUG_ER */
652
653 Value one;
654 value_init(one);
655 value_set_si(one, 1);
656 int len = P->NbRays-1;
657 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
658 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
659 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
660 for (int j = 0; j < P->NbRays; ++j) {
661 if (j == r)
662 continue;
663 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
664 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
665 }
666
667 P = Rays2Polyhedron(M, options->MaxRays);
668 Matrix_Free(M);
669 evalue *EP = barvinok_enumerate_e_with_options(P, exist, nparam, options);
670 Polyhedron_Free(P);
671 value_clear(one);
672
673 return EP;
674 }
675
676 static evalue* enumerate_redundant_ray(Polyhedron *P,
677 unsigned exist, unsigned nparam, barvinok_options *options)
678 {
679 assert(P->NbBid == 0);
680 int nvar = P->Dimension - exist - nparam;
681 Value m;
682 value_init(m);
683
684 for (int r = 0; r < P->NbRays; ++r) {
685 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
686 continue;
687 int i1 = single_param_pos(P, exist, nparam, r);
688 if (i1 == -1)
689 continue;
690 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
691 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
692 continue;
693 int i2 = single_param_pos(P, exist, nparam, r2);
694 if (i2 == -1)
695 continue;
696 if (i1 != i2)
697 continue;
698
699 value_division(m, P->Ray[r][1+nvar+exist+i1],
700 P->Ray[r2][1+nvar+exist+i1]);
701 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
702 /* r2 divides r => r redundant */
703 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
704 value_clear(m);
705 return enumerate_remove_ray(P, r, exist, nparam, options);
706 }
707
708 value_division(m, P->Ray[r2][1+nvar+exist+i1],
709 P->Ray[r][1+nvar+exist+i1]);
710 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
711 /* r divides r2 => r2 redundant */
712 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
713 value_clear(m);
714 return enumerate_remove_ray(P, r2, exist, nparam, options);
715 }
716 }
717 }
718 value_clear(m);
719 return 0;
720 }
721
722 static Polyhedron *upper_bound(Polyhedron *P,
723 int pos, Value *max, Polyhedron **R)
724 {
725 Value v;
726 int r;
727 value_init(v);
728
729 *R = 0;
730 Polyhedron *N;
731 Polyhedron *B = 0;
732 for (Polyhedron *Q = P; Q; Q = N) {
733 N = Q->next;
734 for (r = 0; r < P->NbRays; ++r) {
735 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
736 value_pos_p(P->Ray[r][1+pos]))
737 break;
738 }
739 if (r < P->NbRays) {
740 Q->next = *R;
741 *R = Q;
742 continue;
743 } else {
744 Q->next = B;
745 B = Q;
746 }
747 for (r = 0; r < P->NbRays; ++r) {
748 if (value_zero_p(P->Ray[r][P->Dimension+1]))
749 continue;
750 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
751 if ((!Q->next && r == 0) || value_gt(v, *max))
752 value_assign(*max, v);
753 }
754 }
755 value_clear(v);
756 return B;
757 }
758
759 static evalue* enumerate_ray(Polyhedron *P,
760 unsigned exist, unsigned nparam, barvinok_options *options)
761 {
762 assert(P->NbBid == 0);
763 int nvar = P->Dimension - exist - nparam;
764
765 int r;
766 for (r = 0; r < P->NbRays; ++r)
767 if (value_zero_p(P->Ray[r][P->Dimension+1]))
768 break;
769 if (r >= P->NbRays)
770 return 0;
771
772 int r2;
773 for (r2 = r+1; r2 < P->NbRays; ++r2)
774 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
775 break;
776 if (r2 < P->NbRays) {
777 if (nvar > 0)
778 return enumerate_sum(P, exist, nparam, options);
779 }
780
781 #ifdef DEBUG_ER
782 fprintf(stderr, "\nER: Ray\n");
783 #endif /* DEBUG_ER */
784
785 Value m;
786 Value one;
787 value_init(m);
788 value_init(one);
789 value_set_si(one, 1);
790 int i = single_param_pos(P, exist, nparam, r);
791 assert(i != -1); // for now;
792
793 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
794 for (int j = 0; j < P->NbRays; ++j) {
795 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
796 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
797 }
798 Polyhedron *S = Rays2Polyhedron(M, options->MaxRays);
799 Matrix_Free(M);
800 Polyhedron *D = DomainDifference(P, S, options->MaxRays);
801 Polyhedron_Free(S);
802 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
803 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
804 Polyhedron *R;
805 D = upper_bound(D, nvar+exist+i, &m, &R);
806 assert(D);
807 Domain_Free(D);
808
809 M = Matrix_Alloc(2, P->Dimension+2);
810 value_set_si(M->p[0][0], 1);
811 value_set_si(M->p[1][0], 1);
812 value_set_si(M->p[0][1+nvar+exist+i], -1);
813 value_set_si(M->p[1][1+nvar+exist+i], 1);
814 value_assign(M->p[0][1+P->Dimension], m);
815 value_oppose(M->p[1][1+P->Dimension], m);
816 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
817 P->Ray[r][1+nvar+exist+i]);
818 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
819 // Matrix_Print(stderr, P_VALUE_FMT, M);
820 D = AddConstraints(M->p[0], 2, P, options->MaxRays);
821 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
822 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
823 P->Ray[r][1+nvar+exist+i]);
824 // Matrix_Print(stderr, P_VALUE_FMT, M);
825 S = AddConstraints(M->p[0], 1, P, options->MaxRays);
826 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
827 Matrix_Free(M);
828
829 evalue *EP = barvinok_enumerate_e_with_options(D, exist, nparam, options);
830 Polyhedron_Free(D);
831 value_clear(one);
832 value_clear(m);
833
834 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
835 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, options->MaxRays);
836 else {
837 M = Matrix_Alloc(1, nparam+2);
838 value_set_si(M->p[0][0], 1);
839 value_set_si(M->p[0][1+i], 1);
840 enumerate_vd_add_ray(EP, M, options->MaxRays);
841 Matrix_Free(M);
842 }
843
844 if (!emptyQ(S)) {
845 evalue *E = barvinok_enumerate_e_with_options(S, exist, nparam, options);
846 eadd(E, EP);
847 evalue_free(E);
848 }
849 Polyhedron_Free(S);
850
851 if (R) {
852 assert(nvar == 0);
853 evalue *ER = enumerate_or(R, exist, nparam, options);
854 eor(ER, EP);
855 free_evalue_refs(ER);
856 free(ER);
857 }
858
859 return EP;
860 }
861
862 static evalue* enumerate_vd(Polyhedron **PA,
863 unsigned exist, unsigned nparam, barvinok_options *options)
864 {
865 Polyhedron *P = *PA;
866 int nvar = P->Dimension - exist - nparam;
867 Param_Polyhedron *PP = NULL;
868 Polyhedron *C = Universe_Polyhedron(nparam);
869 Polyhedron *CEq;
870 Matrix *CT;
871 Polyhedron *PR = P;
872 PP = Polyhedron2Param_Polyhedron(PR, C, options);
873 Polyhedron_Free(C);
874
875 int nd;
876 Param_Domain *D, *last;
877 Value c;
878 value_init(c);
879 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
880 ;
881
882 Polyhedron **VD = new Polyhedron *[nd];
883 Polyhedron *TC = true_context(P, C, options->MaxRays);
884 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
885 VD[nd++] = rVD;
886 last = D;
887 END_FORALL_REDUCED_DOMAIN
888 Polyhedron_Free(TC);
889
890 evalue *EP = 0;
891
892 if (nd == 0)
893 EP = evalue_zero();
894
895 /* This doesn't seem to have any effect */
896 if (nd == 1) {
897 Polyhedron *CA = align_context(VD[0], P->Dimension, options->MaxRays);
898 Polyhedron *O = P;
899 P = DomainIntersection(P, CA, options->MaxRays);
900 if (O != *PA)
901 Polyhedron_Free(O);
902 Polyhedron_Free(CA);
903 if (emptyQ(P))
904 EP = evalue_zero();
905 }
906
907 if (PR != *PA)
908 Polyhedron_Free(PR);
909 PR = 0;
910
911 if (!EP && nd > 1) {
912 #ifdef DEBUG_ER
913 fprintf(stderr, "\nER: VD\n");
914 #endif /* DEBUG_ER */
915 for (int i = 0; i < nd; ++i) {
916 Polyhedron *CA = align_context(VD[i], P->Dimension, options->MaxRays);
917 Polyhedron *I = DomainIntersection(P, CA, options->MaxRays);
918
919 if (i == 0)
920 EP = barvinok_enumerate_e_with_options(I, exist, nparam, options);
921 else {
922 evalue *E = barvinok_enumerate_e_with_options(I, exist, nparam,
923 options);
924 eadd(E, EP);
925 evalue_free(E);
926 }
927 Polyhedron_Free(I);
928 Polyhedron_Free(CA);
929 }
930 }
931
932 for (int i = 0; i < nd; ++i)
933 Polyhedron_Free(VD[i]);
934 delete [] VD;
935 value_clear(c);
936
937 if (!EP && nvar == 0) {
938 Value f;
939 value_init(f);
940 Param_Vertices *V, *V2;
941 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
942
943 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
944 bool found = false;
945 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
946 if (V == V2) {
947 found = true;
948 continue;
949 }
950 if (!found)
951 continue;
952 for (int i = 0; i < exist; ++i) {
953 value_oppose(f, V->Vertex->p[i][nparam+1]);
954 Vector_Combine(V->Vertex->p[i],
955 V2->Vertex->p[i],
956 M->p[0] + 1 + nvar + exist,
957 V2->Vertex->p[i][nparam+1],
958 f,
959 nparam+1);
960 int j;
961 for (j = 0; j < nparam; ++j)
962 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
963 break;
964 if (j >= nparam)
965 continue;
966 ConstraintSimplify(M->p[0], M->p[0],
967 P->Dimension+2, &f);
968 value_set_si(M->p[0][0], 0);
969 Polyhedron *para = AddConstraints(M->p[0], 1, P,
970 options->MaxRays);
971 POL_ENSURE_VERTICES(para);
972 if (emptyQ(para)) {
973 Polyhedron_Free(para);
974 continue;
975 }
976 Polyhedron *pos, *neg;
977 value_set_si(M->p[0][0], 1);
978 value_decrement(M->p[0][P->Dimension+1],
979 M->p[0][P->Dimension+1]);
980 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
981 value_set_si(f, -1);
982 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
983 P->Dimension+1);
984 value_decrement(M->p[0][P->Dimension+1],
985 M->p[0][P->Dimension+1]);
986 value_decrement(M->p[0][P->Dimension+1],
987 M->p[0][P->Dimension+1]);
988 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
989 POL_ENSURE_VERTICES(neg);
990 POL_ENSURE_VERTICES(pos);
991 if (emptyQ(neg) && emptyQ(pos)) {
992 Polyhedron_Free(para);
993 Polyhedron_Free(pos);
994 Polyhedron_Free(neg);
995 continue;
996 }
997 #ifdef DEBUG_ER
998 fprintf(stderr, "\nER: Order\n");
999 #endif /* DEBUG_ER */
1000 EP = barvinok_enumerate_e_with_options(para, exist, nparam,
1001 options);
1002 evalue *E;
1003 if (!emptyQ(pos)) {
1004 E = barvinok_enumerate_e_with_options(pos, exist, nparam,
1005 options);
1006 eadd(E, EP);
1007 evalue_free(E);
1008 }
1009 if (!emptyQ(neg)) {
1010 E = barvinok_enumerate_e_with_options(neg, exist, nparam,
1011 options);
1012 eadd(E, EP);
1013 evalue_free(E);
1014 }
1015 Polyhedron_Free(para);
1016 Polyhedron_Free(pos);
1017 Polyhedron_Free(neg);
1018 break;
1019 }
1020 if (EP)
1021 break;
1022 } END_FORALL_PVertex_in_ParamPolyhedron;
1023 if (EP)
1024 break;
1025 } END_FORALL_PVertex_in_ParamPolyhedron;
1026
1027 if (!EP) {
1028 /* Search for vertex coordinate to split on */
1029 /* First look for one independent of the parameters */
1030 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
1031 for (int i = 0; i < exist; ++i) {
1032 int j;
1033 for (j = 0; j < nparam; ++j)
1034 if (value_notzero_p(V->Vertex->p[i][j]))
1035 break;
1036 if (j < nparam)
1037 continue;
1038 value_set_si(M->p[0][0], 1);
1039 Vector_Set(M->p[0]+1, 0, nvar+exist);
1040 Vector_Copy(V->Vertex->p[i],
1041 M->p[0] + 1 + nvar + exist, nparam+1);
1042 value_oppose(M->p[0][1+nvar+i],
1043 V->Vertex->p[i][nparam+1]);
1044
1045 Polyhedron *pos, *neg;
1046 value_set_si(M->p[0][0], 1);
1047 value_decrement(M->p[0][P->Dimension+1],
1048 M->p[0][P->Dimension+1]);
1049 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
1050 value_set_si(f, -1);
1051 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
1052 P->Dimension+1);
1053 value_decrement(M->p[0][P->Dimension+1],
1054 M->p[0][P->Dimension+1]);
1055 value_decrement(M->p[0][P->Dimension+1],
1056 M->p[0][P->Dimension+1]);
1057 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1058 POL_ENSURE_VERTICES(neg);
1059 POL_ENSURE_VERTICES(pos);
1060 if (emptyQ(neg) || emptyQ(pos)) {
1061 Polyhedron_Free(pos);
1062 Polyhedron_Free(neg);
1063 continue;
1064 }
1065 Polyhedron_Free(pos);
1066 value_increment(M->p[0][P->Dimension+1],
1067 M->p[0][P->Dimension+1]);
1068 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1069 #ifdef DEBUG_ER
1070 fprintf(stderr, "\nER: Vertex\n");
1071 #endif /* DEBUG_ER */
1072 pos->next = neg;
1073 EP = enumerate_or(pos, exist, nparam, options);
1074 break;
1075 }
1076 if (EP)
1077 break;
1078 } END_FORALL_PVertex_in_ParamPolyhedron;
1079 }
1080
1081 if (!EP) {
1082 /* Search for vertex coordinate to split on */
1083 /* Now look for one that depends on the parameters */
1084 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
1085 for (int i = 0; i < exist; ++i) {
1086 value_set_si(M->p[0][0], 1);
1087 Vector_Set(M->p[0]+1, 0, nvar+exist);
1088 Vector_Copy(V->Vertex->p[i],
1089 M->p[0] + 1 + nvar + exist, nparam+1);
1090 value_oppose(M->p[0][1+nvar+i],
1091 V->Vertex->p[i][nparam+1]);
1092
1093 Polyhedron *pos, *neg;
1094 value_set_si(M->p[0][0], 1);
1095 value_decrement(M->p[0][P->Dimension+1],
1096 M->p[0][P->Dimension+1]);
1097 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
1098 value_set_si(f, -1);
1099 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
1100 P->Dimension+1);
1101 value_decrement(M->p[0][P->Dimension+1],
1102 M->p[0][P->Dimension+1]);
1103 value_decrement(M->p[0][P->Dimension+1],
1104 M->p[0][P->Dimension+1]);
1105 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1106 POL_ENSURE_VERTICES(neg);
1107 POL_ENSURE_VERTICES(pos);
1108 if (emptyQ(neg) || emptyQ(pos)) {
1109 Polyhedron_Free(pos);
1110 Polyhedron_Free(neg);
1111 continue;
1112 }
1113 Polyhedron_Free(pos);
1114 value_increment(M->p[0][P->Dimension+1],
1115 M->p[0][P->Dimension+1]);
1116 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
1117 #ifdef DEBUG_ER
1118 fprintf(stderr, "\nER: ParamVertex\n");
1119 #endif /* DEBUG_ER */
1120 pos->next = neg;
1121 EP = enumerate_or(pos, exist, nparam, options);
1122 break;
1123 }
1124 if (EP)
1125 break;
1126 } END_FORALL_PVertex_in_ParamPolyhedron;
1127 }
1128
1129 Matrix_Free(M);
1130 value_clear(f);
1131 }
1132
1133 if (CEq)
1134 Polyhedron_Free(CEq);
1135 if (CT)
1136 Matrix_Free(CT);
1137 if (PP)
1138 Param_Polyhedron_Free(PP);
1139 *PA = P;
1140
1141 return EP;
1142 }
1143
1144 evalue* barvinok_enumerate_pip(Polyhedron *P, unsigned exist, unsigned nparam,
1145 unsigned MaxRays)
1146 {
1147 evalue *E;
1148 barvinok_options *options = barvinok_options_new_with_defaults();
1149 options->MaxRays = MaxRays;
1150 E = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1151 barvinok_options_free(options);
1152 return E;
1153 }
1154
1155 evalue *barvinok_enumerate_pip_with_options(Polyhedron *P,
1156 unsigned exist, unsigned nparam, struct barvinok_options *options)
1157 {
1158 int nvar = P->Dimension - exist - nparam;
1159 evalue *EP = evalue_zero();
1160 Polyhedron *Q, *N;
1161
1162 #ifdef DEBUG_ER
1163 fprintf(stderr, "\nER: PIP\n");
1164 #endif /* DEBUG_ER */
1165
1166 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
1167 for (Q = D; Q; Q = N) {
1168 N = Q->next;
1169 Q->next = 0;
1170 evalue *E;
1171 exist = Q->Dimension - nvar - nparam;
1172 E = barvinok_enumerate_e_with_options(Q, exist, nparam, options);
1173 Polyhedron_Free(Q);
1174 eadd(E, EP);
1175 evalue_free(E);
1176 }
1177
1178 return EP;
1179 }
1180
1181 static bool is_single(Value *row, int pos, int len)
1182 {
1183 return First_Non_Zero(row, pos) == -1 &&
1184 First_Non_Zero(row+pos+1, len-pos-1) == -1;
1185 }
1186
1187 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1188 unsigned exist, unsigned nparam, barvinok_options *options);
1189
1190 #ifdef DEBUG_ER
1191 static int er_level = 0;
1192
1193 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1194 unsigned exist, unsigned nparam, barvinok_options *options)
1195 {
1196 fprintf(stderr, "\nER: level %i\n", er_level);
1197
1198 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
1199 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
1200 ++er_level;
1201 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1202 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1203 Polyhedron_Free(P);
1204 --er_level;
1205 return EP;
1206 }
1207 #else
1208 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
1209 unsigned exist, unsigned nparam, barvinok_options *options)
1210 {
1211 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
1212 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
1213 Polyhedron_Free(P);
1214 return EP;
1215 }
1216 #endif
1217
1218 evalue* barvinok_enumerate_e(Polyhedron *P, unsigned exist, unsigned nparam,
1219 unsigned MaxRays)
1220 {
1221 evalue *E;
1222 barvinok_options *options = barvinok_options_new_with_defaults();
1223 options->MaxRays = MaxRays;
1224 E = barvinok_enumerate_e_with_options(P, exist, nparam, options);
1225 barvinok_options_free(options);
1226 return E;
1227 }
1228
1229 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
1230 unsigned exist, unsigned nparam, barvinok_options *options)
1231 {
1232 if (exist == 0) {
1233 Polyhedron *U = Universe_Polyhedron(nparam);
1234 evalue *EP = barvinok_enumerate_with_options(P, U, options);
1235 //char *param_name[] = {"P", "Q", "R", "S", "T" };
1236 //print_evalue(stdout, EP, param_name);
1237 Polyhedron_Free(U);
1238 return EP;
1239 }
1240
1241 int nvar = P->Dimension - exist - nparam;
1242 int len = P->Dimension + 2;
1243
1244 /* for now */
1245 POL_ENSURE_FACETS(P);
1246 POL_ENSURE_VERTICES(P);
1247
1248 if (emptyQ(P))
1249 return evalue_zero();
1250
1251 if (nvar == 0 && nparam == 0) {
1252 evalue *EP = evalue_zero();
1253 barvinok_count_with_options(P, &EP->x.n, options);
1254 if (value_pos_p(EP->x.n))
1255 value_set_si(EP->x.n, 1);
1256 return EP;
1257 }
1258
1259 int r;
1260 for (r = 0; r < P->NbRays; ++r)
1261 if (value_zero_p(P->Ray[r][0]) ||
1262 value_zero_p(P->Ray[r][P->Dimension+1])) {
1263 int i;
1264 for (i = 0; i < nvar; ++i)
1265 if (value_notzero_p(P->Ray[r][i+1]))
1266 break;
1267 if (i >= nvar)
1268 continue;
1269 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
1270 if (value_notzero_p(P->Ray[r][i+1]))
1271 break;
1272 if (i >= nvar + exist + nparam)
1273 break;
1274 }
1275 if (r < P->NbRays) {
1276 evalue *EP = evalue_zero();
1277 value_set_si(EP->x.n, -1);
1278 return EP;
1279 }
1280
1281 int first;
1282 for (r = 0; r < P->NbEq; ++r)
1283 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
1284 break;
1285 if (r < P->NbEq) {
1286 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
1287 exist-first-1) != -1) {
1288 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1289 #ifdef DEBUG_ER
1290 fprintf(stderr, "\nER: Equality\n");
1291 #endif /* DEBUG_ER */
1292 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1293 options);
1294 Polyhedron_Free(T);
1295 return EP;
1296 } else {
1297 #ifdef DEBUG_ER
1298 fprintf(stderr, "\nER: Fixed\n");
1299 #endif /* DEBUG_ER */
1300 if (first == 0)
1301 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1302 options);
1303 else {
1304 Polyhedron *T = Polyhedron_Copy(P);
1305 SwapColumns(T, nvar+1, nvar+1+first);
1306 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1307 options);
1308 Polyhedron_Free(T);
1309 return EP;
1310 }
1311 }
1312 }
1313
1314 Vector *row = Vector_Alloc(len);
1315 value_set_si(row->p[0], 1);
1316
1317 Value f;
1318 value_init(f);
1319
1320 enum constraint* info = new constraint[exist];
1321 for (int i = 0; i < exist; ++i) {
1322 info[i] = ALL_POS;
1323 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1324 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1325 continue;
1326 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
1327 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1328 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1329 continue;
1330 bool lu_parallel = l_parallel ||
1331 is_single(P->Constraint[u]+nvar+1, i, exist);
1332 value_oppose(f, P->Constraint[u][nvar+i+1]);
1333 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
1334 f, P->Constraint[l][nvar+i+1], len-1);
1335 if (!(info[i] & INDEPENDENT)) {
1336 int j;
1337 for (j = 0; j < exist; ++j)
1338 if (j != i && value_notzero_p(row->p[nvar+j+1]))
1339 break;
1340 if (j == exist) {
1341 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
1342 info[i] = (constraint)(info[i] | INDEPENDENT);
1343 }
1344 }
1345 if (info[i] & ALL_POS) {
1346 value_addto(row->p[len-1], row->p[len-1],
1347 P->Constraint[l][nvar+i+1]);
1348 value_addto(row->p[len-1], row->p[len-1], f);
1349 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1350 value_subtract(row->p[len-1], row->p[len-1], f);
1351 value_decrement(row->p[len-1], row->p[len-1]);
1352 ConstraintSimplify(row->p, row->p, len, &f);
1353 value_set_si(f, -1);
1354 Vector_Scale(row->p+1, row->p+1, f, len-1);
1355 value_decrement(row->p[len-1], row->p[len-1]);
1356 Polyhedron *T = AddConstraints(row->p, 1, P, options->MaxRays);
1357 POL_ENSURE_VERTICES(T);
1358 if (!emptyQ(T)) {
1359 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
1360 info[i] = (constraint)(info[i] ^ ALL_POS);
1361 }
1362 //puts("pos remainder");
1363 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1364 Polyhedron_Free(T);
1365 }
1366 if (!(info[i] & ONE_NEG)) {
1367 if (lu_parallel) {
1368 negative_test_constraint(P->Constraint[l],
1369 P->Constraint[u],
1370 row->p, nvar+i, len, &f);
1371 oppose_constraint(row->p, len, &f);
1372 Polyhedron *T = AddConstraints(row->p, 1, P,
1373 options->MaxRays);
1374 POL_ENSURE_VERTICES(T);
1375 if (emptyQ(T)) {
1376 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
1377 info[i] = (constraint)(info[i] | ONE_NEG);
1378 }
1379 //puts("neg remainder");
1380 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1381 Polyhedron_Free(T);
1382 } else if (!(info[i] & ROT_NEG)) {
1383 if (parallel_constraints(P->Constraint[l],
1384 P->Constraint[u],
1385 row->p, nvar, exist)) {
1386 negative_test_constraint7(P->Constraint[l],
1387 P->Constraint[u],
1388 row->p, nvar, exist,
1389 len, &f);
1390 oppose_constraint(row->p, len, &f);
1391 Polyhedron *T = AddConstraints(row->p, 1, P,
1392 options->MaxRays);
1393 POL_ENSURE_VERTICES(T);
1394 if (emptyQ(T)) {
1395 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
1396 info[i] = (constraint)(info[i] | ROT_NEG);
1397 r = l;
1398 }
1399 //puts("neg remainder");
1400 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1401 Polyhedron_Free(T);
1402 }
1403 }
1404 }
1405 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
1406 goto next;
1407 }
1408 }
1409 if (info[i] & ALL_POS)
1410 break;
1411 next:
1412 ;
1413 }
1414
1415 /*
1416 for (int i = 0; i < exist; ++i)
1417 printf("%i: %i\n", i, info[i]);
1418 */
1419 for (int i = 0; i < exist; ++i)
1420 if (info[i] & ALL_POS) {
1421 #ifdef DEBUG_ER
1422 fprintf(stderr, "\nER: Positive\n");
1423 #endif /* DEBUG_ER */
1424 // Eliminate
1425 // Maybe we should chew off some of the fat here
1426 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
1427 for (int j = 0; j < P->Dimension; ++j)
1428 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
1429 Polyhedron *T = Polyhedron_Image(P, M, options->MaxRays);
1430 Matrix_Free(M);
1431 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1432 options);
1433 Polyhedron_Free(T);
1434 value_clear(f);
1435 Vector_Free(row);
1436 delete [] info;
1437 return EP;
1438 }
1439 for (int i = 0; i < exist; ++i)
1440 if (info[i] & ONE_NEG) {
1441 #ifdef DEBUG_ER
1442 fprintf(stderr, "\nER: Negative\n");
1443 #endif /* DEBUG_ER */
1444 Vector_Free(row);
1445 value_clear(f);
1446 delete [] info;
1447 if (i == 0)
1448 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
1449 options);
1450 else {
1451 Polyhedron *T = Polyhedron_Copy(P);
1452 SwapColumns(T, nvar+1, nvar+1+i);
1453 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1454 options);
1455 Polyhedron_Free(T);
1456 return EP;
1457 }
1458 }
1459 for (int i = 0; i < exist; ++i)
1460 if (info[i] & ROT_NEG) {
1461 #ifdef DEBUG_ER
1462 fprintf(stderr, "\nER: Rotate\n");
1463 #endif /* DEBUG_ER */
1464 Vector_Free(row);
1465 value_clear(f);
1466 delete [] info;
1467 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
1468 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
1469 options);
1470 Polyhedron_Free(T);
1471 return EP;
1472 }
1473 for (int i = 0; i < exist; ++i)
1474 if (info[i] & INDEPENDENT) {
1475 Polyhedron *pos, *neg;
1476
1477 /* Find constraint again and split off negative part */
1478
1479 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1480 row, f, true, &pos, &neg)) {
1481 #ifdef DEBUG_ER
1482 fprintf(stderr, "\nER: Split\n");
1483 #endif /* DEBUG_ER */
1484
1485 evalue *EP =
1486 barvinok_enumerate_e_with_options(neg, exist-1, nparam, options);
1487 evalue *E =
1488 barvinok_enumerate_e_with_options(pos, exist, nparam, options);
1489 eadd(E, EP);
1490 evalue_free(E);
1491 Polyhedron_Free(neg);
1492 Polyhedron_Free(pos);
1493 value_clear(f);
1494 Vector_Free(row);
1495 delete [] info;
1496 return EP;
1497 }
1498 }
1499 delete [] info;
1500
1501 Polyhedron *O = P;
1502 Polyhedron *F;
1503
1504 evalue *EP;
1505
1506 EP = enumerate_line(P, exist, nparam, options);
1507 if (EP)
1508 goto out;
1509
1510 EP = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
1511 if (EP)
1512 goto out;
1513
1514 EP = enumerate_redundant_ray(P, exist, nparam, options);
1515 if (EP)
1516 goto out;
1517
1518 EP = enumerate_sure(P, exist, nparam, options);
1519 if (EP)
1520 goto out;
1521
1522 EP = enumerate_ray(P, exist, nparam, options);
1523 if (EP)
1524 goto out;
1525
1526 EP = enumerate_sure2(P, exist, nparam, options);
1527 if (EP)
1528 goto out;
1529
1530 F = unfringe(P, options->MaxRays);
1531 if (!PolyhedronIncludes(F, P)) {
1532 #ifdef DEBUG_ER
1533 fprintf(stderr, "\nER: Fringed\n");
1534 #endif /* DEBUG_ER */
1535 EP = barvinok_enumerate_e_with_options(F, exist, nparam, options);
1536 Polyhedron_Free(F);
1537 goto out;
1538 }
1539 Polyhedron_Free(F);
1540
1541 if (nparam)
1542 EP = enumerate_vd(&P, exist, nparam, options);
1543 if (EP)
1544 goto out2;
1545
1546 if (nvar != 0) {
1547 EP = enumerate_sum(P, exist, nparam, options);
1548 goto out2;
1549 }
1550
1551 assert(nvar == 0);
1552
1553 int i;
1554 Polyhedron *pos, *neg;
1555 for (i = 0; i < exist; ++i)
1556 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
1557 row, f, false, &pos, &neg))
1558 break;
1559
1560 assert (i < exist);
1561
1562 pos->next = neg;
1563 EP = enumerate_or(pos, exist, nparam, options);
1564
1565 out2:
1566 if (O != P)
1567 Polyhedron_Free(P);
1568
1569 out:
1570 value_clear(f);
1571 Vector_Free(row);
1572 return EP;
1573 }
lattice point enumeration library