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