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