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