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