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